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