001 /*
002 * Java Genetic Algorithm Library (jenetics-4.3.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 and space complexity of {@code O(n^2},
186 * where {@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 d[i][j] = dominance.compare(set.get(i), set.get(j));
210 d[j][i] = -d[i][j];
211 }
212 }
213
214 // Compute for each element p the element q that it dominates and the
215 // number of times it is dominated. Using the names as defined in the
216 // referenced paper.
217 final int[] nq = new int[set.size()];
218 final List<IntList> fronts = new ArrayList<>();
219 IntList Fi = new IntList();
220
221 for (int p = 0; p < set.size(); ++p) {
222 final IntList Sp = new IntList();
223 int np = 0;
224
225 for (int q = 0; q < set.size(); ++q) {
226 if (p != q) {
227 // If p dominates q, add q to the set of solutions
228 // dominated by p.
229 if (d[p][q] > 0) {
230 Sp.add(q);
231
232 // Increment the domination counter of p.
233 } else if (d[q][p] > 0) {
234 np += 1;
235 }
236 }
237 }
238
239 // p belongs to the first front.
240 if (np == 0) {
241 Fi.add(p);
242 }
243
244 fronts.add(Sp);
245 nq[p] = np;
246 }
247
248 // Initialize the front counter.
249 int i = 0;
250 final int[] ranks = new int[set.size()];
251 while (!Fi.isEmpty()) {
252 // Used to store the members of the next front.
253 final IntList Q = new IntList();
254
255 for (int p = 0; p < Fi.size(); ++p) {
256 final int fi = Fi.get(p);
257 ranks[fi] = i;
258
259 // Update the dominated counts as compute next front.
260 for (int k = 0, n = fronts.get(fi).size(); k < n; ++k) {
261 final int q = fronts.get(fi).get(k);
262 nq[q] -= 1;
263
264 // q belongs to the next front.
265 if (nq[q] == 0) {
266 Q.add(q);
267 }
268 }
269 }
270
271 ++i;
272 Fi = Q;
273 }
274
275 return ranks;
276 }
277
278 /* *************************************************************************
279 * 'front'
280 * ************************************************************************/
281
282 /**
283 * Return the elements, from the given input {@code set}, which are part of
284 * the pareto front. The {@link Comparable} interface defines the dominance
285 * measure of the elements, used for calculating the pareto front.
286 * <p>
287 * <b>Reference:</b><em>
288 * E. Zitzler and L. Thiele.
289 * Multiobjective Evolutionary Algorithms: A Comparative Case Study
290 * and the Strength Pareto Approach,
291 * IEEE Transactions on Evolutionary Computation, vol. 3, no. 4,
292 * pp. 257-271, 1999.</em>
293 *
294 * @param set the input set
295 * @param <T> the element type
296 * @return the elements which are part of the pareto set
297 * @throws NullPointerException if one of the arguments is {@code null}
298 */
299 public static <T> ISeq<Vec<T>> front(final Seq<? extends Vec<T>> set) {
300 return front(set, Vec::dominance);
301 }
302
303 /**
304 * Return the elements, from the given input {@code set}, which are part of
305 * the pareto front.
306 * <p>
307 * <b>Reference:</b><em>
308 * E. Zitzler and L. Thiele.
309 * Multiobjective Evolutionary Algorithms: A Comparative Case Study
310 * and the Strength Pareto Approach,
311 * IEEE Transactions on Evolutionary Computation, vol. 3, no. 4,
312 * pp. 257-271, 1999.</em>
313 *
314 * @param set the input set
315 * @param dominance the dominance comparator used
316 * @param <T> the element type
317 * @return the elements which are part of the pareto set
318 * @throws NullPointerException if one of the arguments is {@code null}
319 */
320 public static <T> ISeq<T> front(
321 final Seq<? extends T> set,
322 final Comparator<? super T> dominance
323 ) {
324 final MSeq<T> front = MSeq.of(set);
325
326 int n = front.size();
327 int i = 0;
328 while (i < n) {
329 int j = i + 1;
330 while (j < n) {
331 if (dominance.compare(front.get(i), front.get(j)) > 0) {
332 --n;
333 front.swap(j, n);
334 } else if (dominance.compare(front.get(j), front.get(i)) > 0) {
335 --n;
336 front.swap(i, n);
337 --i;
338 break;
339 } else {
340 ++j;
341 }
342 }
343 ++i;
344 }
345
346 return front.subSeq(0, n).copy().toISeq();
347 }
348
349
350 /* *************************************************************************
351 * Common 'dominance' methods.
352 * ************************************************************************/
353
354 /**
355 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
356 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
357 *
358 * @see Vec#dominance(Comparable[], Comparable[])
359 *
360 * @param u the first vector
361 * @param v the second vector
362 * @param <C> the element type of vector <b>u</b> and <b>v</b>
363 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
364 * <b>u</b> and {@code 0} otherwise
365 * @throws NullPointerException if one of the arguments is {@code null}
366 * @throws IllegalArgumentException if {@code u.length != v.length}
367 */
368 public static <C extends Comparable<? super C>> int
369 dominance(final C[] u, final C[] v) {
370 return dominance(u, v, Comparator.naturalOrder());
371 }
372
373 /**
374 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
375 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
376 *
377 * @see Vec#dominance(Object[], Object[], Comparator)
378 *
379 * @param u the first vector
380 * @param v the second vector
381 * @param comparator the element comparator which is used for calculating
382 * the dominance
383 * @param <T> the element type of vector <b>u</b> and <b>v</b>
384 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
385 * <b>u</b> and {@code 0} otherwise
386 * @throws NullPointerException if one of the arguments is {@code null}
387 * @throws IllegalArgumentException if {@code u.length != v.length}
388 */
389 public static <T> int
390 dominance(final T[] u, final T[] v, final Comparator<? super T> comparator) {
391 requireNonNull(comparator);
392 checkLength(u.length, v.length);
393
394 return dominance(
395 u, v, u.length, (a, b, i) -> comparator.compare(a[i], b[i])
396 );
397 }
398
399 /**
400 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
401 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
402 *
403 * @see Vec#dominance(int[], int[])
404 *
405 * @param u the first vector
406 * @param v the second vector
407 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
408 * <b>u</b> and {@code 0} otherwise
409 * @throws NullPointerException if one of the arguments is {@code null}
410 * @throws IllegalArgumentException if {@code u.length != v.length}
411 */
412 public static int dominance(final int[] u, final int[] v) {
413 checkLength(u.length, v.length);
414
415 return dominance(
416 u, v, u.length, (a, b, i) -> Integer.compare(a[i], b[i])
417 );
418 }
419
420 /**
421 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
422 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
423 *
424 * @see Vec#dominance(long[], long[])
425 *
426 * @param u the first vector
427 * @param v the second vector
428 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
429 * <b>u</b> and {@code 0} otherwise
430 * @throws NullPointerException if one of the arguments is {@code null}
431 * @throws IllegalArgumentException if {@code u.length != v.length}
432 */
433 public static int dominance(final long[] u, final long[] v) {
434 checkLength(u.length, v.length);
435
436 return dominance(
437 u, v, u.length, (a, b, i) -> Long.compare(a[i], b[i])
438 );
439 }
440
441 /**
442 * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
443 * <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
444 *
445 * @see Vec#dominance(double[], double[])
446 *
447 * @param u the first vector
448 * @param v the second vector
449 * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
450 * <b>u</b> and {@code 0} otherwise
451 * @throws NullPointerException if one of the arguments is {@code null}
452 * @throws IllegalArgumentException if {@code u.length != v.length}
453 */
454 public static int dominance(final double[] u, final double[] v) {
455 checkLength(u.length, v.length);
456
457 return dominance(
458 u, v, u.length, (a, b, i) -> Double.compare(a[i], b[i])
459 );
460 }
461
462 private static void checkLength(final int i, final int j) {
463 if (i != j) {
464 throw new IllegalArgumentException(format(
465 "Length are not equals: %d != %d.", i, j
466 ));
467 }
468 }
469
470 private static <T> int dominance(
471 final T u,
472 final T v,
473 final int dimension,
474 final ElementComparator<? super T> comparator
475 ) {
476 boolean udominated = false;
477 boolean vdominated = false;
478
479 for (int i = 0; i < dimension; ++i) {
480 final int cmp = comparator.compare(u, v, i);
481
482 if (cmp > 0) {
483 udominated = true;
484 if (vdominated) {
485 return 0;
486 }
487 } else if (cmp < 0) {
488 vdominated = true;
489 if (udominated) {
490 return 0;
491 }
492 }
493 }
494
495 if (udominated == vdominated) {
496 return 0;
497 } else if (udominated) {
498 return 1;
499 } else {
500 return -1;
501 }
502 }
503
504 }
|