001/*
002 * Java Genetic Algorithm Library (jenetics-7.2.0).
003 * Copyright (c) 2007-2023 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 */
020package io.jenetics.ext.moea;
021
022import static java.lang.Double.POSITIVE_INFINITY;
023import static java.lang.String.format;
024import static java.util.Objects.requireNonNull;
025
026import java.util.ArrayList;
027import java.util.Arrays;
028import java.util.Comparator;
029import java.util.List;
030import java.util.function.ToIntFunction;
031
032import io.jenetics.util.BaseSeq;
033import io.jenetics.util.ISeq;
034import io.jenetics.util.MSeq;
035import io.jenetics.util.ProxySorter;
036
037import io.jenetics.ext.internal.util.IntList;
038
039/**
040 * Low-level utility methods for doing pareto-optimal calculations. These methods
041 * are mostly for users who want 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 */
047public 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(BaseSeq, 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[]
076        crowdingDistance(final BaseSeq<? 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(BaseSeq)
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 BaseSeq<? 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.length()];
118                if (set.length() < 3) {
119                        Arrays.fill(result, POSITIVE_INFINITY);
120                } else {
121                        for (int m = 0, d = dimension.applyAsInt(set.get(0)); m < d; ++m) {
122                                final int[] idx = ProxySorter.sort(
123                                        set,
124                                        comparator.ofIndex(m).reversed()
125                                );
126
127                                result[idx[0]] = POSITIVE_INFINITY;
128                                result[idx[set.length() - 1]] = POSITIVE_INFINITY;
129
130                                final T max = set.get(idx[0]);
131                                final T min = set.get(idx[set.length() - 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.length() - 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 BaseSeq<? 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 BaseSeq<? extends T> set,
203                final Comparator<? super T> dominance
204        ) {
205                // Pre-compute the dominance relations.
206                final int[][] d = new int[set.length()][set.length()];
207                for (int i = 0; i < set.length(); ++i) {
208                        for (int j = i + 1; j < set.length(); ++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.length()];
218                final List<IntList> fronts = new ArrayList<>();
219                IntList Fi = new IntList();
220
221                for (int p = 0; p < set.length(); ++p) {
222                        final IntList Sp = new IntList();
223                        int np = 0;
224
225                        for (int q = 0; q < set.length(); ++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.length()];
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                                // The updates 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 BaseSeq<? 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 BaseSeq<? 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        /**
471         * Calculates the <a href="https://en.wikipedia.org/wiki/Pareto_efficiency">
472         *     <b>Pareto Dominance</b></a> of the two vectors <b>u</b> and <b>v</b>.
473         *
474         * @since 5.2
475         *
476         * @param u the first vector
477         * @param v the second vector
478         * @param dimensions the number of vector elements
479         * @param comparator the comparator used for comparing the vector elements
480         * @param <V> the vector type
481         * @return {@code 1} if <b>u</b> ≻ <b>v</b>, {@code -1} if <b>v</b> ≻
482         *         <b>u</b> and {@code 0} otherwise
483         * @throws NullPointerException if one of the arguments is {@code null}
484         */
485        public static <V> int dominance(
486                final V u,
487                final V v,
488                final int dimensions,
489                final ElementComparator<? super V> comparator
490        ) {
491                boolean udominated = false;
492                boolean vdominated = false;
493
494                for (int i = 0; i < dimensions; ++i) {
495                        final int cmp = comparator.compare(u, v, i);
496
497                        if (cmp > 0) {
498                                if (vdominated) {
499                                        return 0;
500                                } else {
501                                        udominated = true;
502                                }
503                        } else if (cmp < 0) {
504                                if (udominated) {
505                                        return 0;
506                                } else {
507                                        vdominated = true;
508                                }
509                        }
510                }
511
512                if (udominated == vdominated) {
513                        return 0;
514                } else if (udominated) {
515                        return 1;
516                } else {
517                        return -1;
518                }
519        }
520
521}