[pre-commit.ci] pre-commit autoupdate

.pre-commit-config.yaml

@@ -7,7 +7,7 @@ ci:
 
 repos:
   - repo: https://github.com/lorenzwalthert/precommit
-    rev: v0.4.3.9020
+    rev: v0.4.3.9021
     hooks:
       - id: parsable-R
       - id: no-browser-statement
@@ -50,14 +50,14 @@ repos:
       - id: tox-ini-fmt
 
   - repo: https://github.com/tox-dev/pyproject-fmt
-    rev: v2.16.1
+    rev: v2.18.1
     hooks:
       - id: pyproject-fmt
         additional_dependencies: ["tox>=4.12.1"]
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
     # Ruff version.
-    rev: v0.15.1
+    rev: v0.15.6
     hooks:
       # Run the formatter.
       - id: ruff-format

TODO

@@ -18,11 +18,7 @@ Computation, pp. 465–472 (2012)]
        and https://dl.acm.org/doi/pdf/10.1145/3205455.3205469 and http://www.cmap.polytechnique.fr/~nikolaus.hansen/proceedings/2016/GECCO/proceedings/p613.pdf
  - Online nondominated filter
    - [ ] [A Fast Incremental BSP Tree Archive for Non-dominated Points](https://doi.org/10.1007/978-3-319-54157-0_18)
- - [ ] Faster nondominated algorithm (is_nondominated, filter_dominated) for d>3:
-    - Kung: https://doi.org/10.1145/321906.321910
-    - [Bentley: Multidimensional Divide-and-Conquer](https://dl.acm.org/doi/pdf/10.1145/358841.358850)
-    - C++ implementation: https://github.com/TLDart/nondLib/blob/main/nondlib.hpp <https://dl.acm.org/doi/10.1145/3449726.3462737>
-    - Faster for very many points: [Worst-case I/O-efficient Skyline Algorithms](https://www.cse.cuhk.edu.hk/~taoyf/paper/tods12-sky.pdf)
+ - [ ] Faster nondominated algorithm (is_nondominated, filter_dominated) for very many points: [Worst-case I/O-efficient Skyline Algorithms](https://www.cse.cuhk.edu.hk/~taoyf/paper/tods12-sky.pdf)
  - [ ] Faster nondominated sorting: https://github.com/mbuzdalov/non-dominated-sorting       https://github.com/KernelA/nds-py
  <https://discourse.julialang.org/t/fast-non-dominated-sorting/24164/20>
  <https://ctlab.itmo.ru/~mbuzdalov/ec-2018/05-nds.pdf>
@@ -50,3 +46,7 @@ Computation, pp. 465–472 (2012)]
  - [X] move eaftools/eaf code here.
  - [x] move hypervolume/trunk code here.
  - [x] move hypervolume testsuite here.
+ - [x] Faster nondominated algorithm (is_nondominated, filter_dominated) for d>3:
+    - Kung et al.: https://doi.org/10.1145/321906.321910
+    - [Bentley: Multidimensional Divide-and-Conquer](https://dl.acm.org/doi/pdf/10.1145/358841.358850)
+    - [C++ implementation](https://github.com/TLDart/nondLib/blob/main/nondlib.hpp) and [paper about the implementation](https://dl.acm.org/doi/10.1145/3449726.3462737)

c/Makefile

@@ -106,7 +106,9 @@ HDRS =  avl.h                                                                \
         libmoocore-config.h	                                             \
         mt19937/mt19937.h                                                    \
         nondominated.h                                                       \
-        pow_int.h                                                            \
+        nondominated_kung.h                                                  \
+	nth_element.h                                                        \
+	pow_int.h                                                            \
         rng.h                                                                \
         sort.h                                                               \
         timer.h                                                              \
@@ -202,24 +204,27 @@ $(EXE_FILES): cmdline.o io.o
 
 avl.o: avl.h
 cmdline.o: cmdline.h io.h
-dominatedsets.o: cmdline.h io.h nondominated.h sort.h avl_tiny.h epsilon.h
+dominatedsets.o: cmdline.h io.h nondominated.h epsilon.h
 eaf.o: eaf.h io.h bit_array.h cvector.h
 eaf3d.o: eaf.h io.h bit_array.h cvector.h avl.h
 eaf_main.o: cmdline.h io.h eaf.h bit_array.h cvector.h
-epsilon.o: cmdline.h io.h epsilon.h nondominated.h sort.h avl_tiny.h
+epsilon.o: cmdline.h io.h nondominated.h epsilon.h
 hv.o: hv.h hv_priv.h hvc4d_priv.h sort.h libmoocore-config.h
 hv3dplus.o: hv_priv.h sort.h hv3d_priv.h avl_tiny.h
 hv4d.o: hv4d_priv.h hv_priv.h sort.h
-hv_contrib.o: hv.h libmoocore-config.h nondominated.h sort.h avl_tiny.h
+hv_contrib.o: hv.h libmoocore-config.h nondominated.h
 hvapprox.o: hvapprox.h pow_int.h rng.h mt19937/mt19937.h
 hvc3d.o: hv_priv.h sort.h hv3d_priv.h avl_tiny.h
-igd.o: cmdline.h io.h igd.h nondominated.h sort.h avl_tiny.h pow_int.h
+igd.o: cmdline.h io.h nondominated.h igd.h pow_int.h
 io.o: io.h io_priv.h
-main-hv.o: cmdline.h io.h hv.h timer.h libmoocore-config.h nondominated.h sort.h avl_tiny.h hvapprox.h
+main-hv.o: cmdline.h io.h hv.h timer.h libmoocore-config.h nondominated.h hvapprox.h
 main-hvapprox.o: cmdline.h io.h hvapprox.h
-ndsort.o: cmdline.h io.h nondominated.h sort.h avl_tiny.h hv.h
-nondominated.o : cmdline.h io.h nondominated.h avl_tiny.h
-pareto.o: nondominated.h sort.h avl_tiny.h
+ndsort.o: cmdline.h io.h nondominated.h hv.h
+nondominated.h: avl_tiny.h nondominated_kung.h
+nondominated.o : cmdline.h io.h nondominated.h
+nondominated_kung.h: nth_element.h
+nth_element.h: sort.h
+pareto.o: nondominated.h
 rng.o: rng.h mt19937/mt19937.h ziggurat_constants.h
 timer.o: timer.h
 $(OBJS): common.h gcc_attribs.h config.h maxminclamp.h

c/NEWS.md

@@ -4,6 +4,7 @@
 
  * `fpli_hv()` calculates 4D contributions as the base case, which is
    significantly faster for more than four dimensions. (Andreia P. Guerreiro)
+ * nondominated_kung.h: New.
 
 ## 0.18
 

c/nondominated.h

@@ -346,14 +346,16 @@ find_nondominated_set_2d_(const double * restrict points, size_t size,
    A different implementation is available from Duarte M. Dias, Alexandre
    D. Jesus, Luís Paquete, A software library for archiving nondominated
    points, GECCO 2021. https://github.com/TLDart/nondLib/blob/main/nondlib.hpp
+
+   rows: must be sorted by cmp_ppdouble_asc_rev_3d().
+   max_dim is only used if the rows pointers have been shifted.
 */
 static inline size_t
-find_nondominated_set_3d_helper(const double * restrict points, size_t size,
+find_nondominated_set_3d_helper(const double * restrict points, const double ** rows,
+                                size_t size, dimension_t max_dim,
                                 const bool keep_weakly, bool * restrict nondom)
 {
-    ASSUME(size >= 2);
-    const double ** rows = generate_row_pointers_asc_rev_3d(points, size);
-
+    ASSUME(size > 1);
     avl_tree_t tree;
     avl_init_tree(&tree, qsort_cmp_pdouble_asc_x_nonzero);
     avl_node_t * tnodes = malloc((size+1) * sizeof(*tnodes));
@@ -370,9 +372,9 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
     const double * last_dom = NULL;
     const double * restrict pk = rows[0];
     do {
-        const double * restrict pj = rows[j];
-        DEBUG3(printf_point("pj: ", pj, 3, "\n"));
         bool dominated;
+        const double * restrict pj = rows[j];
+        DEBUG2(printf_point("pj = [ ", pj, 3, " ], "));
         if (pk[0] > pj[0] || pk[1] > pj[1]) {
             avl_node_t * nodeaux;
             int res = avl_search_closest(&tree, pj, &nodeaux);
@@ -382,14 +384,14 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
                 assert(prev[0] != sentinel[0]);
                 assert(prev[0] <= pj[0]);
                 dominated = prev[1] <= pj[1];
-                DEBUG3(printf_point("res > 0: prev: ", prev, 3, "\n"));
+                DEBUG2(printf_point("res > 0: prev: ", prev, 3, "\n"));
                 nodeaux = nodeaux->next;
             } else if (nodeaux->prev) { // nodeaux goes after pj, so move to the next one.
                 const double * restrict prev = nodeaux->prev->item;
                 assert(prev[0] != sentinel[0]);
                 assert(prev[0] <= pj[0]);
                 dominated = prev[1] <= pj[1];
-                DEBUG3(printf_point("res < 0: prev: ", prev, 3, "\n"));
+                DEBUG2(printf_point("res < 0: prev: ", prev, 3, "\n"));
             } else {
                 dominated = false;
             }
@@ -399,7 +401,7 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
                 assert(pj[0] <= point[0]);
                 // Delete everything in the tree that is dominated by pj.
                 while (pj[1] <= point[1]) {
-                    DEBUG3(printf_point("delete point: ", point, 3, "\n"));
+                    DEBUG2(printf_point("delete point: ", point, 3, "\n"));
                     assert(pj[0] <= point[0]);
                     nodeaux = nodeaux->next;
                     point = nodeaux->item;
@@ -408,7 +410,7 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
                        rebalance. */
                     avl_unlink_node(&tree, nodeaux->prev);
                 }
-                DEBUG3((point == sentinel)
+                DEBUG2((point == sentinel)
                        ? printf_point("insert before sentinel: ", sentinel, 2, "\n")
                        : printf_point("insert before point: ", point, 3, "\n"));
                 (++node)->item = pj;
@@ -427,12 +429,12 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
                     // or pk was dominated, then this one is also dominated.
                     || last_dom == pk;
             }
-            DEBUG3(printf_point("weakly dominated by pk: ", pk, 3, "\n"));
+            DEBUG2(printf_point("weakly dominated by pk: ", pk, 3, "\n"));
         }
         if (dominated) { // pj is dominated by a point in the tree or by prev.
-            /* Map the order in p[], which is sorted, to the original order in
-               points. */
-            size_t pos_last_dom = row_index_from_ptr(points, pj, 3);
+            /* Map the order in rows[], which is sorted, to the original order
+               in points. */
+            size_t pos_last_dom = row_index_from_shifted_ptr(points, pj, 3, max_dim);
             if (unlikely(nondom == NULL)) {
                 // In this context, it means "position of the first dominated solution found".
                 new_size = pos_last_dom;
@@ -449,25 +451,35 @@ find_nondominated_set_3d_helper(const double * restrict points, size_t size,
 
 early_end:
     free(tnodes);
-    free(rows);
     return new_size;
 }
 
 static inline size_t
 find_dominated_3d_(const double * restrict points, size_t size, bool keep_weakly)
 {
-    return find_nondominated_set_3d_helper(points, size, keep_weakly, /* nondom=*/NULL);
+    const double ** rows = generate_row_pointers_asc_rev_3d(points, size);
+    size_t pos = keep_weakly
+        ? find_nondominated_set_3d_helper(points, rows, size, 3, true, /* nondom=*/NULL)
+        : find_nondominated_set_3d_helper(points, rows, size, 3, false, /* nondom=*/NULL);
+    free(rows);
+    return pos;
 }
 
 /*
    Store which points are nondominated in nondom and return the number of
    nondominated points.
 */
 static inline size_t
-find_nondominated_set_3d_(const double * restrict points, size_t size, const bool keep_weakly, bool * restrict nondom)
+find_nondominated_set_3d_(const double * restrict points, size_t size,
+                          const bool keep_weakly, bool * restrict nondom)
 {
+    const double ** rows = generate_row_pointers_asc_rev_3d(points, size);
     ASSUME(nondom != NULL);
-    return find_nondominated_set_3d_helper(points, size, keep_weakly, nondom);
+    size_t new_size = keep_weakly
+        ? find_nondominated_set_3d_helper(points, rows, size, 3, true, nondom)
+        : find_nondominated_set_3d_helper(points, rows, size, 3, false, nondom);
+    free(rows);
+    return new_size;
 }
 
 /**
@@ -601,6 +613,8 @@ find_dominated_point_(const double * restrict points, size_t size, dimension_t d
     }
 }
 
+#include "nondominated_kung.h"
+
 /**
    Store which points are nondominated in nondom and return the number of
    nondominated points.
@@ -618,14 +632,19 @@ find_nondominated_set_(const double * restrict points, size_t size, dimension_t
     ASSUME(minmax != NULL);
     ASSUME(nondom != NULL);
 
-    if (dim <= 3) {
+    /* FIXME: In theory, for large dim and small n, brute-force may be faster
+       than Kung's algorithm. How to set this threshold requires further
+       investigation.  */
+    if (dim <= 3 || size > KUNG_SMALL_THRESHOLD) {
         const double * pp = force_agree_minimize(points, size, &dim, agree, minmax);
         ASSUME(dim >= 2);
         size_t res;
         if (dim == 2) {
             res = find_nondominated_set_2d_(pp, size, keep_weakly, nondom);
-        } else {
+        } else if (dim == 3) {
             res = find_nondominated_set_3d_(pp, size, keep_weakly, nondom);
+        } else {
+            res = find_nondominated_set_agree_kung(pp, size, dim, keep_weakly, nondom);
         }
         if (pp != points)
             free((void *) pp);
@@ -654,6 +673,7 @@ find_nondominated_set_(const double * restrict points, size_t size, dimension_t
           unreachable();
     }
 }
+#undef KUNG_SMALL_THRESHOLD
 
 static inline size_t
 find_dominated_point_agree(const double * restrict points, size_t size, dimension_t dim,