diff --git a/README.md b/README.md index 423fa2e7d..532511fec 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,8 @@ Mickaël Binois, Leonardo C.T. Bezerra, Fergus Rooney, - [Lennart Schäpermeier](https://schaepermeier.github.io). + [Lennart Schäpermeier](https://schaepermeier.github.io), + Jonathan Fieldsend. --------------------------------------- diff --git a/TODO-archiving b/TODO-archiving new file mode 100644 index 000000000..27a0f29de --- /dev/null +++ b/TODO-archiving @@ -0,0 +1,14 @@ + +# Unbounded Archive + - [ ] .add_many() is currently a loop in Python and it is not available in C. + First, it could loop in C. + Also, it could filter before adding unless the user specifies "nondominated=True". + If it is nondominated, it could do a faster merge operation. + + - Benchmark against: https://github.com/alandefreitas/pareto + + +# Bounded size archives. + + - [ ] Hypervolume archiver + - [ ] MGA diff --git a/bibkeys.txt b/bibkeys.txt index 0a0ff5285..e9e8edff5 100644 --- a/bibkeys.txt +++ b/bibkeys.txt @@ -1,3 +1,4 @@ +AraSei1989rst AugBadBroZit2009gecco BeuFonLopPaqVah09:tec BezLopStu2017assessment @@ -27,6 +28,7 @@ HanJas1998 HerWer1987tabucol IshHeSha2019regular IshMasTanNoj2015igd +JasLus2018ndtree Jen03 JohAraMcGSch1991 KnoCor2002cec diff --git a/c/Makefile b/c/Makefile index 81922b8d6..4f2780728 100644 --- a/c/Makefile +++ b/c/Makefile @@ -72,6 +72,7 @@ SRCS = avl.c \ igd.c \ io.c \ libutil.c \ + main-archiver.o \ main-hv.c \ main-hvapprox.c \ main-nondominated.c \ @@ -81,11 +82,12 @@ SRCS = avl.c \ pareto.c \ rng.c \ timer.c \ + treap2d_fixed_size.c \ + treap_archive.c \ whv.c \ whv_hype.c HDRS = avl.h \ - avl_tiny.h \ bit_array.h \ cmdline.h \ common.h \ @@ -112,6 +114,9 @@ HDRS = avl.h \ rng.h \ sort.h \ timer.h \ + treap2d.h \ + treap2d_dlnode.h \ + treap2d_fixed_size.h \ whv.h \ whv_hype.h \ ziggurat_constants.h @@ -160,7 +165,7 @@ EXE_CFLAGS += $(SANITIZERS) $(OPT_CFLAGS) $(MARCH_FLAGS) $(WARN_CFLAGS) \ # List of test target base names -TEST_NAMES := eaf nondominated hv hvapprox igd epsilon dominatedsets ndsort +TEST_NAMES := archiver eaf nondominated hv hvapprox igd epsilon dominatedsets ndsort ASAN_TARGETS := $(addprefix asan-,$(TEST_NAMES)) UBSAN_TARGETS:= $(addprefix ubsan-,$(TEST_NAMES)) TEST_TARGETS := $(addprefix test-,$(TEST_NAMES)) @@ -192,14 +197,15 @@ $(ALL_NAMES): %: $(BINDIR)/%$(EXE) include libhv.mk -$(BINDIR)/dominatedsets$(EXE): dominatedsets.o nondominated.o +$(BINDIR)/archiver$(EXE): main-archiver.o io.o treap_archive.o nondominated.o treap2d_fixed_size.o +$(BINDIR)/dominatedsets$(EXE): dominatedsets.o nondominated.o treap2d_fixed_size.o $(BINDIR)/eaf$(EXE): eaf.o eaf_main.o eaf3d.o avl.o -$(BINDIR)/epsilon$(EXE): epsilon.o nondominated.o +$(BINDIR)/epsilon$(EXE): epsilon.o nondominated.o treap2d_fixed_size.o $(BINDIR)/hv$(EXE): main-hv.o timer.o $(LIBHV_OBJS) -$(BINDIR)/hvapprox$(EXE): main-hvapprox.o timer.o hvapprox.o rng.o mt19937/mt19937.o nondominated.o -$(BINDIR)/igd$(EXE): igd.o nondominated.o -$(BINDIR)/ndsort$(EXE): ndsort.o pareto.o nondominated.o $(LIBHV_OBJS) -$(BINDIR)/nondominated$(EXE): main-nondominated.o timer.o nondominated.o +$(BINDIR)/hvapprox$(EXE): main-hvapprox.o timer.o hvapprox.o rng.o mt19937/mt19937.o nondominated.o treap2d_fixed_size.o +$(BINDIR)/igd$(EXE): igd.o nondominated.o treap2d_fixed_size.o +$(BINDIR)/ndsort$(EXE): ndsort.o pareto.o $(LIBHV_OBJS) +$(BINDIR)/nondominated$(EXE): main-nondominated.o timer.o nondominated.o treap2d_fixed_size.o $(EXE_FILES): cmdline.o io.o $(call MKDIR, $(BINDIR)/) @@ -214,22 +220,25 @@ 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 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 +hv3dplus.o: hv_priv.h sort.h hv3d_priv.h treap2d.h hv4d.o: hv4d_priv.h hv_priv.h sort.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 +hvc3d.o: hv_priv.h sort.h hv3d_priv.h treap2d.h igd.o: cmdline.h io.h nondominated.h igd.h pow_int.h io.o: io.h io_priv.h +main-archiver.o: cmdline.h io.h treap_archive.h treap2d.h nd_tree.h libmoocore-config.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 main-nondominated.o : cmdline.h io.h nondominated.h nondominated_kung.h ndsort.o: cmdline.h io.h nondominated.h hv.h -nondominated.h: avl_tiny.h nondominated_kung.h -nondominated.o: nondominated.h nondominated_kung.h -pareto.o: nondominated.h +nondominated.h: nondominated_kung.h +nondominated.o: nondominated.h nondominated_kung.h treap2d_fixed_size.h +pareto.o: nondominated.h treap2d_fixed_size.h rng.o: rng.h mt19937/mt19937.h ziggurat_constants.h timer.o: timer.h +treap2d_fixed_size.o: treap2d_fixed_size.h +treap2d_fixed_size.h: treap2d.h $(OBJS): common.h gcc_attribs.h config.h maxminclamp.h shlibs: clean $(SHLIB) $(CXXSHLIB) diff --git a/c/NEWS.md b/c/NEWS.md index 317b10c1d..c1aa478ff 100644 --- a/c/NEWS.md +++ b/c/NEWS.md @@ -6,6 +6,9 @@ * main-nondominated.c: Renamed from nondominated.c. * hv and hvapprox executables now use -U instead of -u for --union to be consistent with nondominated. + * treap2d.h: New. + * pareto.c, nondominated.c, hv3d_priv.h: Use Treap instead of AVL-tree. + * avl_tiny.h: Delete. ## 0.19.2 diff --git a/c/README.md b/c/README.md index 4dd10db26..ff31252e8 100644 --- a/c/README.md +++ b/c/README.md @@ -26,6 +26,7 @@ The [**moocore** project](https://github.com/multi-objective/moocore/) collects * Computation of the Empirical Attainment Function. The empirical attainment function (EAF) describes the probabilistic distribution of the outcomes obtained by a stochastic algorithm in the objective space. + * Archivers that allow efficiently storing a nondominated set from points given one-by-one as a sequence. The goal of **moocore** is to provide an implementation of these functions that is efficient, thoroughly tested, well-documented, multi-platform (Windows, Linux, MacOS), multi-language (C, R, and Python, so far), with minimal dependencies and reusable from other libraries and packages, thus providing an easy-to-use and high-performance building block. diff --git a/c/archive.h b/c/archive.h new file mode 100644 index 000000000..e79ea3eec --- /dev/null +++ b/c/archive.h @@ -0,0 +1,110 @@ +#ifndef _ARCHIVE_H_ +#define _ARCHIVE_H_ + +#include "treap_archive.h" +#include "nd_tree.h" + +enum unbounded_archive_type_t { + UNBOUNDED_ARCHIVE_TREAP, + UNBOUNDED_ARCHIVE_NDTREE, +}; + +typedef struct UnboundedArchive { + enum unbounded_archive_type_t archive_type; + union { + TreapArchive treap; + NDTreeArchive ndtree; + }; +} UnboundedArchive; + +#define UNBOUNDED_ARCHIVE_DISPATCH(SELF, TREAP_CODE, NDTREE_CODE) \ + switch ((SELF)->archive_type) { \ + case UNBOUNDED_ARCHIVE_TREAP: do { TREAP_CODE; } while(0); break; \ + case UNBOUNDED_ARCHIVE_NDTREE: do { NDTREE_CODE; } while(0); break; \ + default: unreachable(); \ + } + + +static inline UnboundedArchive * +archive_new(enum unbounded_archive_type_t type) +{ + UnboundedArchive * arch = malloc(sizeof(*arch)); + arch->archive_type = type; + return arch; + +} +static inline UnboundedArchive * +archive_ndtree_new(dimension_t dim, NDTreeConfig config) +{ + UnboundedArchive * archive = archive_new(UNBOUNDED_ARCHIVE_NDTREE); + ndtree_init(&archive->ndtree, dim, config.max_children, config.max_bucket_size, config.allow_duplicates); + return archive; +} + +static inline UnboundedArchive * +archive_treap_new(void) +{ + UnboundedArchive * archive = archive_new(UNBOUNDED_ARCHIVE_TREAP); + treap_archive_init(&archive->treap); + return archive; +} + +/** + +*/ +static inline int +archive_add(UnboundedArchive * restrict arch, const double * restrict z, const void * restrict x) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + return treap_archive_add_discard(&arch->treap, z, x), + return ndtree_add_discard(&arch->ndtree, z, x)); +} + +static inline size_t +archive_total_size(const UnboundedArchive *arch) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + return treap_archive_total_size(&arch->treap), + return ndtree_total_size(&arch->ndtree)); +} + +static inline size_t +archive_unique_size(const UnboundedArchive *arch) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + return treap_archive_unique_size(&arch->treap), + return ndtree_unique_size(&arch->ndtree)); +} + +static inline void +archive_get_contents(const UnboundedArchive *arch, double * z_list, const void ** x_list) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + treap_archive_get_contents(&arch->treap, z_list, x_list), + ndtree_get_contents(&arch->ndtree, z_list, x_list)); +} + +static inline void +archive_free(UnboundedArchive *arch) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + treap_archive_free_nodes(&arch->treap), + ndtree_free_nodes(&arch->ndtree)); + free(arch); +} + +static inline void +archive_print_stats(FILE * stream, const UnboundedArchive *arch) +{ + UNBOUNDED_ARCHIVE_DISPATCH( + arch, + (void)0, + ndtree_print_stats(stream, &arch->ndtree)); +} + +#endif // _ARCHIVE_H_ diff --git a/c/archiving.h b/c/archiving.h new file mode 100644 index 000000000..224770906 --- /dev/null +++ b/c/archiving.h @@ -0,0 +1,19 @@ +/***************************************************************************** + + Public definitions that are common to all types of archives. + + *****************************************************************************/ +#ifndef _ARCHIVING_H_ +#define _ARCHIVING_H_ + +enum archive_insert_result_t { + ARCHIVE_INSERT_UKNOWN_ERROR = -4, + ARCHIVE_INSERT_X_MUST_BE_NOT_NULL = -3, + ARCHIVE_INSERT_X_MUST_BE_NULL = -2, + ARCHIVE_INSERT_MEMORY_ERROR = -1, + ARCHIVE_INSERT_REJECTED = 0, + ARCHIVE_INSERT_ACCEPTED = 1, + ARCHIVE_INSERT_DUPLICATED = 2, +}; + +#endif // _ARCHIVING_H_ diff --git a/c/archiving_priv.h b/c/archiving_priv.h new file mode 100644 index 000000000..b60af703e --- /dev/null +++ b/c/archiving_priv.h @@ -0,0 +1,26 @@ +#ifndef _ARCHIVING_PRIV_H_ +#define _ARCHIVING_PRIV_H_ + +#include +#include +#include + +#include "archiving.h" + +/** + Returns 0 if OK, negative values if there is a problem. +*/ +static inline int +wrong_x_is_null(signed char * restrict x_is_null, const void * restrict x) +{ + if (*x_is_null == 0) + return x != NULL ? 0 : ARCHIVE_INSERT_X_MUST_BE_NOT_NULL; + if (*x_is_null == 1) + return x == NULL ? 0 : ARCHIVE_INSERT_X_MUST_BE_NULL; + + *x_is_null = (x == NULL); + return 0; +} + + +#endif // _ARCHIVING_PRIV_H_ diff --git a/c/avl_tiny.h b/c/avl_tiny.h deleted file mode 100644 index 94b64d6c5..000000000 --- a/c/avl_tiny.h +++ /dev/null @@ -1,375 +0,0 @@ -#ifndef _AVL_TINY_H_ -#define _AVL_TINY_H_ -/*----------------------------------------------------------------------------- - - The following is a reduced version of the AVL-tree library used here - according to the terms of the LGPL-2.1-or-later. See the copyright notice - below. - -*/ - -/***************************************************************************** - - avl.h - Source code for the AVL-tree library. - - Copyright (C) 1998 Michael H. Buselli - Copyright (C) 2000-2002 Wessel Dankers - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - Augmented AVL-tree. Original by Michael H. Buselli . - - Modified by Wessel Dankers to add a bunch of bloat to - the sourcecode, change the interface and squash a few bugs. - Mail him if you find new bugs. - -*****************************************************************************/ - -/* User supplied function to compare two items like strcmp() does. - * For example: cmp(a,b) will return: - * -1 if a < b - * 0 if a = b - * 1 if a > b - */ -// "const void" here is actually "const avl_item_t" -typedef int (*avl_compare_t)(const void *, const void *); - -typedef struct avl_tree_t { - avl_node_t *head; - avl_node_t *tail; - avl_node_t *top; - avl_compare_t cmp; -} avl_tree_t; - - -/***************************************************************************** - - avl.c - Source code for the AVL-tree library. - -*****************************************************************************/ - -static void avl_rebalance(avl_tree_t *, avl_node_t *); - -#define NODE_DEPTH(n) ((n) ? (n)->depth : 0) -#define L_DEPTH(n) (NODE_DEPTH((n)->left)) -#define R_DEPTH(n) (NODE_DEPTH((n)->right)) -#define CALC_DEPTH(n) (unsigned char)((L_DEPTH(n)>R_DEPTH(n)?L_DEPTH(n):R_DEPTH(n)) + 1) - -static int avl_check_balance(const avl_node_t * restrict avlnode) { - int d = R_DEPTH(avlnode) - L_DEPTH(avlnode); - return d < -1 ? -1 : d > 1 ? 1 : 0; -} - -/* Searches for a node with the key closest (or equal) to the given item. - If avlnode is not NULL, *avlnode will be set to the node found or NULL - if the tree is empty. Return values: - -1 if the returned node is greater - 0 if the returned node is equal or if the tree is empty - 1 if the returned node is smaller - O(lg n) */ -static int -avl_search_closest(const avl_tree_t *avltree, const avl_item_t *item, avl_node_t **avlnode) -{ - avl_node_t *node = avltree->top; - assert(node); - avl_compare_t cmp = avltree->cmp; - while(true) { - int c = cmp(item, node->item); - if(c < 0) { - if(node->left) - node = node->left; - else - return *avlnode = node, -1; - } else if(c > 0) { - if(node->right) - node = node->right; - else - return *avlnode = node, 1; - } else { - return *avlnode = node, 0; - } - } -} - -/* Initializes a new tree for elements that will be ordered using - * the supplied strcmp()-like function. - * Returns the value of avltree (even if it's NULL). - * O(1) */ -static avl_tree_t * -avl_init_tree(avl_tree_t *rc, avl_compare_t cmp) { - rc->head = NULL; - rc->tail = NULL; - rc->top = NULL; - rc->cmp = cmp; - return rc; -} - -/* Reinitializes the tree structure for reuse. Nothing is free()d. - * Compare and freeitem functions are left alone. - * O(1) */ -static inline void -avl_clear_tree(avl_tree_t *avltree) { - avltree->top = avltree->head = avltree->tail = NULL; -} - -static void -avl_clear_node(avl_node_t *newnode) { - newnode->left = newnode->right = NULL; - newnode->depth = 1; -} - -/* Insert a node in an empty tree. If avlnode is NULL, the tree will be - * cleared and ready for re-use. - * If the tree is not empty, the old nodes are left dangling. - * O(1) */ -static avl_node_t * -avl_insert_top(avl_tree_t *avltree, avl_node_t *newnode) { - avl_clear_node(newnode); - newnode->prev = newnode->next = newnode->parent = NULL; - avltree->head = avltree->tail = avltree->top = newnode; - return newnode; -} - -static avl_node_t * -avl_insert_after(avl_tree_t *avltree, avl_node_t *node, avl_node_t *newnode); - -/* Insert a node before another node. Returns the new node. - * If old is NULL, the item is appended to the tree. - * O(lg n) */ -static avl_node_t * -avl_insert_before(avl_tree_t *avltree, avl_node_t *node, avl_node_t *newnode) { - if(!node) - return avltree->tail - ? avl_insert_after(avltree, avltree->tail, newnode) - : avl_insert_top(avltree, newnode); - - if(node->left) - return avl_insert_after(avltree, node->prev, newnode); - - assert (node); - assert (!node->left); - - avl_clear_node(newnode); - - newnode->next = node; - newnode->parent = node; - - newnode->prev = node->prev; - if(node->prev) - node->prev->next = newnode; - else - avltree->head = newnode; - node->prev = newnode; - - node->left = newnode; - avl_rebalance(avltree, node); - return newnode; -} - -/* Insert a node after another node. Returns the new node. - * If old is NULL, the item is prepended to the tree. - * O(lg n) */ -static avl_node_t * -avl_insert_after(avl_tree_t *avltree, avl_node_t *node, avl_node_t *newnode) { - if(!node) - return avltree->head - ? avl_insert_before(avltree, avltree->head, newnode) - : avl_insert_top(avltree, newnode); - - if(node->right) - return avl_insert_before(avltree, node->next, newnode); - - assert (node); - assert (!node->right); - - avl_clear_node(newnode); - - newnode->prev = node; - newnode->parent = node; - - newnode->next = node->next; - if(node->next) - node->next->prev = newnode; - else - avltree->tail = newnode; - node->next = newnode; - - node->right = newnode; - avl_rebalance(avltree, node); - return newnode; -} - -/* Deletes the node from the tree. Does not delete the item at that node. - * The item of the node may be freed before calling avl_unlink_node. - * (In other words, it is not referenced by this function.) - * O(lg n) */ -static void -avl_unlink_node(avl_tree_t *avltree, avl_node_t *avlnode) { - avl_node_t *parent; - avl_node_t **superparent; - avl_node_t *subst, *left, *right; - avl_node_t *balnode; - - if(avlnode->prev) - avlnode->prev->next = avlnode->next; - else - avltree->head = avlnode->next; - - if(avlnode->next) - avlnode->next->prev = avlnode->prev; - else - avltree->tail = avlnode->prev; - - parent = avlnode->parent; - - superparent = parent - ? avlnode == parent->left ? &parent->left : &parent->right - : &avltree->top; - - left = avlnode->left; - right = avlnode->right; - if(!left) { - *superparent = right; - if(right) - right->parent = parent; - balnode = parent; - } else if(!right) { - *superparent = left; - left->parent = parent; - balnode = parent; - } else { - subst = avlnode->prev; - if(subst == left) { - balnode = subst; - } else { - balnode = subst->parent; - balnode->right = subst->left; - if(balnode->right) - balnode->right->parent = balnode; - subst->left = left; - left->parent = subst; - } - subst->right = right; - subst->parent = parent; - right->parent = subst; - *superparent = subst; - } - - avl_rebalance(avltree, balnode); -} - -/* - * avl_rebalance: - * Rebalances the tree if one side becomes too heavy. This function - * assumes that both subtrees are AVL-trees with consistant data. The - * function has the additional side effect of recalculating the count of - * the tree at this node. It should be noted that at the return of this - * function, if a rebalance takes place, the top of this subtree is no - * longer going to be the same node. - */ -static void -avl_rebalance(avl_tree_t *avltree, avl_node_t *avlnode) { - avl_node_t *child; - avl_node_t *gchild; - avl_node_t *parent; - avl_node_t **superparent; - - parent = avlnode; - - while(avlnode) { - parent = avlnode->parent; - - superparent = parent - ? avlnode == parent->left ? &parent->left : &parent->right - : &avltree->top; - - switch(avl_check_balance(avlnode)) { - case -1: - child = avlnode->left; - if(L_DEPTH(child) >= R_DEPTH(child)) { - avlnode->left = child->right; - if(avlnode->left) - avlnode->left->parent = avlnode; - child->right = avlnode; - avlnode->parent = child; - *superparent = child; - child->parent = parent; - avlnode->depth = CALC_DEPTH(avlnode); - child->depth = CALC_DEPTH(child); - } else { - gchild = child->right; - avlnode->left = gchild->right; - if(avlnode->left) - avlnode->left->parent = avlnode; - child->right = gchild->left; - if(child->right) - child->right->parent = child; - gchild->right = avlnode; - if(gchild->right) - gchild->right->parent = gchild; - gchild->left = child; - if(gchild->left) - gchild->left->parent = gchild; - *superparent = gchild; - gchild->parent = parent; - avlnode->depth = CALC_DEPTH(avlnode); - child->depth = CALC_DEPTH(child); - gchild->depth = CALC_DEPTH(gchild); - } - break; - case 1: - child = avlnode->right; - if(R_DEPTH(child) >= L_DEPTH(child)) { - avlnode->right = child->left; - if(avlnode->right) - avlnode->right->parent = avlnode; - child->left = avlnode; - avlnode->parent = child; - *superparent = child; - child->parent = parent; - avlnode->depth = CALC_DEPTH(avlnode); - child->depth = CALC_DEPTH(child); - } else { - gchild = child->left; - avlnode->right = gchild->left; - if(avlnode->right) - avlnode->right->parent = avlnode; - child->left = gchild->right; - if(child->left) - child->left->parent = child; - gchild->left = avlnode; - if(gchild->left) - gchild->left->parent = gchild; - gchild->right = child; - if(gchild->right) - gchild->right->parent = gchild; - *superparent = gchild; - gchild->parent = parent; - avlnode->depth = CALC_DEPTH(avlnode); - child->depth = CALC_DEPTH(child); - gchild->depth = CALC_DEPTH(gchild); - } - break; - default: - avlnode->depth = CALC_DEPTH(avlnode); - } - avlnode = parent; - } -} - -/*----------------------------------------------------------------------------- - end of functions from AVL-tree library. -******************************************************************************/ -#endif diff --git a/c/common.h b/c/common.h index 4e9725bfd..c9d0a0cec 100644 --- a/c/common.h +++ b/c/common.h @@ -142,4 +142,31 @@ euclidean_norm(const double * restrict w, dimension_t dim) return sqrt(norm); } +_attr_pure_func static inline double +squared_distance(const double * restrict a, const double * restrict b, dimension_t dim) +{ + ASSUME(dim >= 2); + double distance = 0; + for (dimension_t i = 0; i < dim; i++) { + double diff = a[i] - b[i]; + distance += diff * diff; + } + ASSUME(distance >= 0); // GCC is not able to infer this. + return distance; +} + +static inline size_t +which_max(const double *v, size_t dim) +{ + double max = v[0]; + size_t pos = 0; + for (size_t i = 1; i < dim; i++) { + if (v[i] > max) { + max = v[i]; + pos = i; + } + } + return pos; +} + #endif /* !MOOCORE_COMMON_H_ */ diff --git a/c/flex_bucket.h b/c/flex_bucket.h new file mode 100644 index 000000000..eccf0ce5d --- /dev/null +++ b/c/flex_bucket.h @@ -0,0 +1,104 @@ +#ifndef FLEX_BUCKET_H_ +#define FLEX_BUCKET_H_ + +#include "solutions_list.h" + +typedef struct FlexBucket { + size_t n_unique; + size_t n_total; + size_t cap; + SolutionsList solutions[]; // List of solutions (each list tracks duplicates). +} FlexBucket; + +static FlexBucket * +flex_bucket_new(size_t cap) +{ + cap = cap > 0 ? cap : 4; + FlexBucket * bucket = malloc(sizeof*bucket + cap * sizeof(SolutionsList)); + if (bucket) { + bucket->cap = cap; + bucket->n_unique = 0; + bucket->n_total = 0; + } + return bucket; +} + +static FlexBucket * +flex_bucket_reserve(FlexBucket * bucket, size_t add_cap) +{ + if (bucket == NULL) + return flex_bucket_new(MAX(add_cap, (size_t)4)); // Use same default cap within flex_bucket_new() + + size_t new_cap = bucket->n_unique + add_cap; + if (new_cap <= bucket->cap) + return bucket; + new_cap = MAX(new_cap, bucket->cap * 2U); + FlexBucket * new_bucket = realloc(bucket, sizeof *bucket + new_cap * sizeof(SolutionsList)); + if (!new_bucket) + return NULL; + new_bucket->cap = new_cap; + return new_bucket; +} + +/** + Copy n_src elements of src to the bucket, possibly re-allocating its memory. + + On return, for each element of src, src.n=0, src[i].x.list = NULL. +*/ +static size_t +flex_bucket_append_move(FlexBucket ** bucket_p, SolutionsList *src, uint32_t n_src) +{ + FlexBucket * bucket = flex_bucket_reserve(*bucket_p, n_src); + if (!bucket) + return 0; + memcpy(&bucket->solutions[bucket->n_unique], src, n_src * sizeof *src); + size_t total = 0; + for (uint32_t i = 0; i < n_src; i++) { + total += src[i].x.n; + // Reset src to detect double copy/free. + src[i].x.n = 0; + src[i].z = NULL; + src[i].x.list = NULL; + } + bucket->n_total += total; + bucket->n_unique += n_src; + *bucket_p = bucket; + return total; +} + +static inline size_t +flex_bucket_unique_len(const FlexBucket * bucket) +{ + return bucket->n_unique; +} + +static inline size_t +flex_bucket_total_len(const FlexBucket * bucket) +{ + return bucket->n_total; +} + +static inline void +flex_bucket_get_unique_vectors(const FlexBucket * bucket, const double ** z_list) +{ + for (size_t i = 0; i < bucket->n_unique; i++, z_list++) + z_list[0] = bucket->solutions[i].z; +} + +static inline void +flex_bucket_get_x_values(const FlexBucket * bucket, const void ** x_list) +{ + for (size_t i = 0; i < bucket->n_unique; i++) + for (size_t j = 0; j < bucket->solutions[i].x.n; j++, x_list++) + x_list[0] = bucket->solutions[i].x.list[j]; +} + +static void +flex_bucket_free(FlexBucket * bucket) +{ + for (size_t i = 0; i < bucket->n_unique; i++) + solutions_list_free(&bucket->solutions[i]); + free(bucket); +} + +#endif // FLEX_BUCKET_H_ diff --git a/c/hv3d_priv.h b/c/hv3d_priv.h index 4168a0779..32f677ccc 100644 --- a/c/hv3d_priv.h +++ b/c/hv3d_priv.h @@ -1,33 +1,7 @@ #ifndef _HV3D_PRIV_H #define _HV3D_PRIV_H -typedef const double avl_item_t; -typedef struct avl_node_t { - struct avl_node_t *next; - struct avl_node_t *prev; - struct avl_node_t *parent; - struct avl_node_t *left; - struct avl_node_t *right; - avl_item_t * item; - dlnode_t * dlnode; - unsigned char depth; -} avl_node_t; - -#include "avl_tiny.h" - -static inline const double * -node_point(const avl_node_t * restrict node) -{ - return node->item; -} - -static inline avl_node_t * -new_avl_node(dlnode_t * restrict p, avl_node_t * restrict node) -{ - node->dlnode = p; - node->item = p->x; - return node; -} +#include "treap2d_dlnode.h" /* Used by hvc3d.c and hv3dplus.c. @@ -61,77 +35,53 @@ hv3d_preprocessing(dlnode_t * restrict list, size_t n) assert(list+1 == list->next[0]); assert(list+2 == list->prev[0]); - avl_tree_t tree; - avl_init_tree(&tree, qsort_cmp_pdouble_asc_y_des_x_nonzero); - avl_node_t * tnodes = malloc((n+2) * sizeof(*tnodes)); - + Treap2D * tree = treap2d_new(n+2); + assert(tree != NULL); // At the top we insert the first point, which is never dominated. dlnode_t * p = (list+1)->next[0]; - avl_node_t * nodeaux = new_avl_node(p, tnodes); - avl_insert_top(&tree, nodeaux); + // Before the top node, we insert list + 1, sentinel 2 (ref[0], -INF), + // After the top node, we insert list, sentinel 1 (-INF, ref[1]) + treap2d_init_with_sentinels(tree, p, list+1, list); set_delimiters(p, list+1, list); - // After the top node, we insert sentinel 1 (-INF, ref[1]) - avl_node_t * node = new_avl_node(list, tnodes + 1); - avl_insert_after(&tree, nodeaux, node); - // Before the top node, we insert sentinel 2 (ref[0], -INF) - node = new_avl_node(list + 1, tnodes + 2); - avl_insert_before(&tree, nodeaux, node); - set_delimiters(p, nodeaux->prev->dlnode, nodeaux->next->dlnode); - + _attr_maybe_unused dlnode_t * prev_p = p; + double pk0 = p->x[0], pk1 = p->x[1], _attr_maybe_unused pk2 = p->x[2]; const dlnode_t * stop = list+2; - p = p->next[0]; - while (p != stop) { - const double * px = p->x; - const double * prev_x; - // == 1 means that nodeaux goes before p, so move to the next one. - if (avl_search_closest(&tree, px, &nodeaux) == 1) { - prev_x = node_point(nodeaux); - nodeaux = nodeaux->next; + for (p = p->next[0]; p != stop; p = p->next[0]) { + const double pj0 = p->x[0], pj1 = p->x[1], pj2 = p->x[2]; + if ((pk0 > pj0) | (pk1 > pj1)) { + TreapNode *pred = treap2d_find_le(tree, pj1); + assert(pred != NULL); + const double * prev_x = treap2d_get_dlnode(pred)->x; + if (prev_x[0] <= pj0) { + // pj is dominated by a point in the tree. +#ifdef HVC_ONLY + if (all_equal_double(prev_x, p->x, 3)) + treap2d_get_dlnode(pred)->ignore = true; // It will have zero hvc. +#endif + remove_from_z(p); + continue; + } + // pj is NOT dominated + TreapNode * prev, *next; + treap2d_insert_and_displace_get_bounds(tree, p, &prev, &next); + // Check that the data structure is properly setup. + assert(treap2d_get_dlnode(prev)->x[0] > pj0 && treap2d_get_dlnode(prev)->x[1] < pj1); + assert(treap2d_get_dlnode(next)->x[0] < pj0 && treap2d_get_dlnode(next)->x[1] > pj1); + set_delimiters(p, treap2d_get_dlnode(prev), treap2d_get_dlnode(next)); + + pk0 = pj0; pk1 = pj1; pk2 = pj2; + prev_p = p; } else { - prev_x = node_point(nodeaux->prev); - } - assert(node_point(nodeaux)[1] > px[1] // node_point(nodeaux) comes after px. - || (node_point(nodeaux)[1] == px[1] && node_point(nodeaux)[0] < px[0])); - assert(prev_x[1] <= px[1]); - assert(prev_x[2] <= px[2]); - if (prev_x[0] <= px[0]) { // px is dominated by a point in the tree. + // pj is dominated by a previous point. #ifdef HVC_ONLY - if (all_equal_double(prev_x, px, 3)) - nodeaux->prev->dlnode->ignore = true; // It will have zero hvc. + if (pk0 == pj0 && pk1 == pj1 && pk2 == pj2) + prev_p->ignore = true; // It will have zero hvc. #endif remove_from_z(p); - } else if (node_point(nodeaux)[1] == px[1]) { // px is dominated by a point in the tree. - // FIXME: If the points were ordered by asc x we would only need the first condition. - assert(node_point(nodeaux)[0] < px[0]); - remove_from_z(p); - } else { - assert(node_point(nodeaux)[1] >= px[1]); - if (prev_x[1] == px[1]){ - nodeaux = nodeaux->prev; - prev_x = node_point(nodeaux->prev); - } - assert(prev_x[1] < px[1]); - // Delete everything in the tree that is dominated by p. - while (node_point(nodeaux)[0] >= px[0]) { - assert(node_point(nodeaux)[1] >= px[1]); - nodeaux = nodeaux->next; - /* FIXME: A possible speed up is to delete without rebalancing - the tree because avl_insert_before() will rebalance, but we - need to know which is the highest node that needs - rebalancing. */ - avl_unlink_node(&tree, nodeaux->prev); - } - node = new_avl_node(p, node + 1); - avl_insert_before(&tree, nodeaux, node); - // Check if the data structure is properly setup - assert(node->prev->dlnode->x[0] > p->x[0] && node->prev->dlnode->x[1] < p->x[1]); - assert(node->next->dlnode->x[0] < p->x[0] && node->next->dlnode->x[1] > p->x[1]); - set_delimiters(p, node->prev->dlnode, node->next->dlnode); } - p = p->next[0]; } - free(tnodes); + treap2d_free(tree); #undef set_delimiters } diff --git a/c/io.c b/c/io.c index 201dd679d..94e6ecba6 100644 --- a/c/io.c +++ b/c/io.c @@ -86,6 +86,7 @@ skip_comment_line(FILE *instream) return 0; } + int fread_double(FILE *instream, double *number) { @@ -163,6 +164,70 @@ fread_int(FILE *instream, int *number) return 1; // Success ! } +/** + This function is used to read one point at a time. + + When *dim == 0, any empty/comment lines are skipped. Then a complete line of + double values is read, stored in point, and *dim is set to the number of + values read, the function returns 1. + + When *dim > 0, if a sequence of empty/comment lines is found, it is skipped and the function returns 0. + Otherwise, a line of *dim double values is read and stored in point, and the function returns 1. + + If EOF is found, return READ_INPUT_FILE_EMPTY. + + If *dim > 0 but a different number of values is found, return READ_INPUT_WRONG_INITIAL_DIM. + + If something different from a number is found, return ERROR_CONVERSION. + + If the number to read is larger than 128 characters, return ERROR_COLUMNS. + +*/ +int +fread_next_double_point(FILE * restrict input, double * restrict point, dimension_t * restrict dim) +{ + assert(input != NULL && point != NULL && dim != NULL); + int status = skip_comment_line(input); + bool found_empty_line = (status == 1); + while (status == 1) + status = skip_comment_line(input); + + if (unlikely(status == EOF)) + return READ_INPUT_FILE_EMPTY; + + dimension_t expected = *dim; + if (expected > 0 && found_empty_line) + return 0; + + // A non-empty line starts here. + // skip_comment_line() has pushed its first character back. + dimension_t n = 0; + do { + if (unlikely(n == MOOCORE_DIMENSION_MAX)) + return ERROR_COLUMNS; + if (unlikely(expected > 0 && n == expected)) + return READ_INPUT_WRONG_INITIAL_DIM; + + double value; + status = fread_double(input, &value); + if (status == 0) + return ERROR_CONVERSION; + if (status == EOF) + break; // Completed line. + + point[n++] = value; + status = skip_comment_line(input); + } while (status == 0); + + if (expected > 0 && n != expected) + return READ_INPUT_WRONG_INITIAL_DIM; + + if (expected == 0) + *dim = n; + + return 1; +} + #define objective_t int #define fread_objective_t fread_int #define read_objective_t_data read_int_data diff --git a/c/io.h b/c/io.h index 5e09a755d..9e3652bfc 100644 --- a/c/io.h +++ b/c/io.h @@ -44,6 +44,9 @@ int write_sets(FILE *outfile, const double *data, int ncols, const int *cumsizes int write_sets_filtered (FILE * outfile, const double * data, int ncols, const int * cumsizes, int nruns, const boolvec * write_p); +int fread_double(FILE *instream, double *number); +int fread_next_double_point(FILE * restrict input, double * restrict point, dimension_t * restrict dim); + static inline const int * read_minmax (const char *str, int *nobj_p) diff --git a/c/libhv.mk b/c/libhv.mk index 7e63a566b..4dcea48d0 100644 --- a/c/libhv.mk +++ b/c/libhv.mk @@ -1,5 +1,5 @@ # -*- Makefile-gmake -*- -LIBHV_SRCS = hv.c hv3dplus.c hv4d.c hvc3d.c hv_contrib.c nondominated.o +LIBHV_SRCS = hv.c hv3dplus.c hv4d.c hvc3d.c hv_contrib.c nondominated.o treap2d_fixed_size.o LIBHV_HDRS = hv.h hv_priv.h hv3d_priv.h hv4d_priv.h hvc4d_priv.h libmoocore-config.h LIBHV_OBJS = $(LIBHV_SRCS:.c=.o) HV_LIB = fpli_hv.a diff --git a/c/main-archiver.c b/c/main-archiver.c new file mode 100644 index 000000000..47220b9c1 --- /dev/null +++ b/c/main-archiver.c @@ -0,0 +1,172 @@ +/************************************************************************* + + archiver: Archive one point at a time. + + --------------------------------------------------------------------- + + Copyright (c) 2026 + Manuel Lopez-Ibanez + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + ---------------------------------------------------------------------- + + Relevant literature: + +*************************************************************************/ +#include "config.h" +#include +#include +#include +#include +#include +#include // for getopt() +#include // for getopt_long() + +#include "io.h" +#include "archive.h" +#define CMDLINE_COPYRIGHT_YEARS "2026" +#define CMDLINE_AUTHORS "Manuel Lopez-Ibanez \n" +#include "cmdline.h" + +static bool union_flag = false; +static bool always_nd_tree_flag = false; +static int verbose_flag = 1; +static char *suffix = NULL; + +static void usage(void) +{ + printf("\n" + "Usage: %s [OPTIONS] [FILE...]\n\n", program_invocation_short_name); + + printf( +"Calculate the hypervolume of each input set of each FILE. \n" +"With no FILE, or when FILE is -, read standard input.\n\n" + +"Options:\n" +OPTION_HELP_STR +OPTION_VERSION_STR +" -v, --verbose print some information (time, maximum, etc). \n" +" --nd-tree Always use ND-Tree for any dimension.\n" +OPTION_UNION_STR +" -s, --suffix=STRING Create an output file for each input file by appending\n" +" this suffix. This is ignored when reading from stdin. \n" +" If missing, output is sent to stdout. \n" +"\n"); +} + +// FIXME: Add options to print quality metrics after archiving each input point. +static bool +read_and_archive(FILE * input) +{ + dimension_t dim = 0; + double point[MOOCORE_DIMENSION_MAX+1]; + // Read the first line as the dimension/template point. + int status = fread_next_double_point(input, point, &dim); + if (status != 1) + fatal_error("reading the first point"); + if (dim < 2) + fatal_error("At least 2 dimensions are required"); + + const NDTreeConfig config = { .max_children = 0, .max_bucket_size = 0, .allow_duplicates = true}; + const bool use_treap = (dim == 2 && !always_nd_tree_flag); + UnboundedArchive * archive = use_treap + ? archive_treap_new() + : archive_ndtree_new(dim, config); + + size_t points_read = 0; + // Read point until EOF or empty line. + do { + if (archive_add(archive, point, /*x=*/NULL) < 0) + fatal_error("Memory error!"); + // FIXME: check that the archive remains nondominated + points_read++; + status = fread_next_double_point(input, point, &dim); + } while (status == 1 || (union_flag && status == 0)); + + if (status != READ_INPUT_FILE_EMPTY && status != 0) + fatal_error("reading point from input"); + + // We completed one data set + size_t total_points = archive_total_size(archive); + // Print result + if (verbose_flag >= 2) { + fprintf(stdout, + "# Points read: %zu\n" + "# Total points: %zu\n" + "# Unique points: %zu\n", + points_read, + total_points, + archive_unique_size(archive)); + } + double * z_list = malloc(sizeof(*z_list) * dim * total_points); + archive_get_contents(archive, z_list, NULL); + // FIXME: Add an option to sort the output. + int n_points = (int) total_points; + // FIXME: write_sets should take size_t * or unsigned *; + write_sets(stdout, z_list, dim, &n_points, 1, /*prefix=*/NULL); + free(z_list); + if (verbose_flag >= 2) + archive_print_stats(stderr, archive); + archive_free(archive); + return (status == 0); +} + + +int main(int argc, char *argv[]) +{ + // See the man page for getopt_long for an explanation of these fields. + static const char short_options[] = "hVvqUs:N"; + static const struct option long_options[] = { + {"help", no_argument, NULL, 'h'}, + {"version", no_argument, NULL, 'V'}, + {"verbose", no_argument, NULL, 'v'}, + {"quiet", no_argument, NULL, 'q'}, + {"union", no_argument, NULL, 'U'}, + {"suffix", required_argument, NULL, 's'}, + {"nd-tree", no_argument, NULL, 'N'}, + {NULL, 0, NULL, 0} /* marks end of list */ + }; + + set_program_invocation_short_name(argv[0]); + + int opt; /* it's actually going to hold a char. */ + int longopt_index; + while (0 < (opt = getopt_long (argc, argv, short_options, + long_options, &longopt_index))) { + switch (opt) { + case 'N': // --nd-tree + always_nd_tree_flag = true; + break; + + case 's': // --suffix + suffix = optarg; + break; + + case 'q': // --quiet + verbose_flag = 0; + break; + + case 'v': // --verbose + verbose_flag = 2; + break; + + case 'U': // --union + union_flag = true; + break; + + default: + default_cmdline_handler(opt); + } + } + + int numfiles = argc - optind; + if (numfiles < 1) { // Read stdin. + while (read_and_archive(stdin)) { + fprintf(stdout, "\n"); + } + } + return EXIT_SUCCESS; +} diff --git a/c/nd_tree.h b/c/nd_tree.h new file mode 100644 index 000000000..369d15288 --- /dev/null +++ b/c/nd_tree.h @@ -0,0 +1,1603 @@ +/****************************************************************************** + + nd_tree.h: ND-Tree Nondominated Archive. + + --------------------------------------------------------------------- + Copyright (C) 2026 + Jonathan Fieldsend + Manuel Lopez-Ibanez + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + --------------------------------------------------------------------- + + Literature: + + [1] Jaszkiewicz, A., Lust, T.: ND-tree-based update: A fast algorithm for the + dynamic nondominance problem. IEEE Transactions on Evolutionary + Computation 22(5), 778-791 (2018). + https://doi.org/10.1109/TEVC.2018.2799684 + + [2] Jaszkiewicz, A., Zielniewicz, P.: Improving the efficiency of the + distance-based hypervolume estimation using ND-Tree. IEEE Transactions on + Evolutionary Computation 29(3), 726-733 (2025). + https://doi.org/10.1109/TEVC.2024.3391857 + + [3] Lang, B.: Space-partitioned ND-Trees for the dynamic nondominance problem. + IEEE Transactions on Evolutionary Computation 26(5), 1004-1014 (2022). + https: //doi.org/10.1109/TEVC.2022.3145631 + + [4] Fieldsend, J.E.: Data structures for non-dominated sets: implementations + and empirical assessment of two decades of advances. In: Coello Coello, + C.A. (ed.) Proceedings of the Genetic and Evolutionary Computation + Conference, GECCO 2020, pp. 489–497. ACM Press, New York, NY (2020). + https://doi.org/10.1145/3377930.3390150 + + [5] Fieldsend, J.E.: On the active use of an ND-Tree-based archive for + multi-objective optimisation. In: Legrand, P., et al. (eds.) EA 2023: + Artificial Evolution, pp. 1-14. Lecture Notes in Computer Science, + Springer (2023). https://doi.org/10.1007/ 978-3-031-42616-2_1 + + + FIXME: The code below does not implement all the improvements described in the + above papers. + + *****************************************************************************/ +#ifndef ND_TREE_ARCHIVE_H +#define ND_TREE_ARCHIVE_H + +#include +#include +#include +#include // memcpy + +#include "sort.h" +#include "maxminclamp.h" +#include "nondominated.h" +#include "archiving_priv.h" +#include "solutions_list.h" +#include "flex_bucket.h" + +/*** + FIXME: The code would be faster if max_bucket_size and/or max_children were + compile-time constants. We could reserve space for 20 max_children but only + use dim+1. +*/ + +/** + ND-tree node. Leaves own buckets; internal nodes own children. The + ideal/nadir/midpoint arrays are conservative bounds used to prune dominance + checks before descending into subtrees. +*/ +typedef struct NDTreeNode { + struct NDTreeNode *parent; + uint8_t n; + uint8_t cap; + // FIXME: Use an anonymous union + SolutionsList *bucket; // List of solutions (each list tracks duplicates). + struct NDTreeNode *children; + // (ideal, nadir) define a hyper-rectangle that contains all points in the + // subtree. The hyper-rectangle is not necessarily tight. + double *ideal; + double *nadir; + double *midpoint; + // This is total_size of subtree (including duplicates). Leaf nodes have + // coverage == 0 and we use node->n. + uint32_t coverage; +} NDTreeNode; + +typedef struct NDTreeConfig { + uint8_t max_children; + uint8_t max_bucket_size; + bool allow_duplicates; +} NDTreeConfig; + +typedef struct NDTreeArchive { + /* It would be faster if this was not a pointer, but then we cannot + relocate the tree in memory, which is needed for NDS-Forest. */ + NDTreeNode * root; + dimension_t dim; + NDTreeConfig config; + signed char x_is_null; + uint64_t revision; // To detect modifications of the tree while iterating. +} NDTreeArchive; + + +static inline SolutionsList * +ndtree_bucket_new(NDTreeNode *node, uint8_t cap) +{ + node->n = 0; + node->cap = cap; + node->bucket = malloc(sizeof(*node->bucket) * node->cap); + return node->bucket; +} + +static inline void +ndtree_bucket_free(NDTreeNode *node) +{ + for (int i = 0; i < node->n; i++) + solutions_list_free(&node->bucket[i]); + free(node->bucket); + node->bucket = NULL; + node->n = 0; +} + +_attr_pure_func +static inline bool +ndtree_node_is_leaf(const NDTreeNode *node) +{ + return node->coverage == 0; +} + +static void +ndtree_node_set_bounds(NDTreeNode * restrict node, const double * restrict z, dimension_t dim) +{ + ASSUME(dim >= 2); + memcpy(node->ideal, z, dim * sizeof(*z)); + memcpy(node->nadir, z, dim * sizeof(*z)); + memcpy(node->midpoint, z, dim * sizeof(*z)); +} + +static inline bool +update_ideal(NDTreeNode * restrict node, const double * restrict z, dimension_t dim) +{ + ASSUME(dim > 1); + bool any_changed = false; + double * restrict ideal = node->ideal; + for (dimension_t i = 0; i < dim; i++) { + if (z[i] < ideal[i]) { + ideal[i] = z[i]; + any_changed = true; + } + } + return any_changed; +} + +static inline bool +update_nadir(NDTreeNode * restrict node, const double * restrict z, dimension_t dim) +{ + ASSUME(dim > 1); + bool any_changed = false; + double * restrict nadir = node->nadir; + for (dimension_t i = 0; i < dim; i++) { + if (z[i] > nadir[i]) { + nadir[i] = z[i]; + any_changed = true; + } + } + return any_changed; +} + +static inline void +update_midpoint(NDTreeNode * restrict node, dimension_t dim) +{ + ASSUME(dim > 1); + const double * restrict ideal = node->ideal; + const double * restrict nadir = node->nadir; + double * restrict midpoint = node->midpoint; + for (dimension_t i = 0; i < dim; i++) { + // With doubles, we do not expect the sum to overflow and this is + // faster than ideal + (nadir - ideal)/2. + midpoint[i] = 0.5 * (ideal[i] + nadir[i]); + } +} + +static inline bool +ndtree_node_update_bounds(NDTreeNode * restrict node, + const double * restrict z, dimension_t dim, + bool * restrict ideal_changed, bool * restrict nadir_changed) +{ + if (*ideal_changed) + *ideal_changed = update_ideal(node, z, dim); + if (*nadir_changed) + *nadir_changed = update_nadir(node, z, dim); + + if (!*ideal_changed && !*nadir_changed) + return false; + update_midpoint(node, dim); + return true; +} + +static void +ndtree_node_update_bounds_up(NDTreeNode * restrict node, const double * restrict z, dimension_t dim) +{ + ASSUME(dim >= 2); + // Iterative instead of recursive because it is faster and helps inlining. + bool ideal_changed = true, nadir_changed = true; + do { + if (!ndtree_node_update_bounds(node, z, dim, &ideal_changed, &nadir_changed)) + return; + node = node->parent; + } while (node); +} + +static void +ndtree_node_propagate_bounds_up(NDTreeNode *node, dimension_t dim) +{ + ASSUME(dim >= 2); + // Propagate the new bounds up the tree. + bool ideal_changed = true, nadir_changed = true; + while (node->parent) { + if (ideal_changed) + ideal_changed = update_ideal(node->parent, node->ideal, dim); + if (nadir_changed) + nadir_changed = update_nadir(node->parent, node->nadir, dim); + if (!ideal_changed && !nadir_changed) + return; + update_midpoint(node->parent, dim); + node = node->parent; + } +} + + +/** + Compute the bounds for a leaf node by looking at the z vectors of its buckets. +*/ +static void +ndtree_node_compute_bounds_of_leaf(NDTreeNode *node, dimension_t dim) +{ + assert(ndtree_node_is_leaf(node)); + assert(node->n > 0); + ndtree_node_set_bounds(node, node->bucket[0].z, dim); + bool any_changed = false; + for (uint8_t i = 1; i < node->n; i++) { + any_changed |= update_ideal(node, node->bucket[i].z, dim); + any_changed |= update_nadir(node, node->bucket[i].z, dim); + } + if (any_changed) + update_midpoint(node, dim); +} + +/** + For a non-leaf node, whose children are all leaves, (re-)compute the bounds + of each child. This function does not propagate the bounds up the tree. + See ndtree_node_propagate_bounds_up() +*/ +static void +ndtree_node_compute_bounds_of_children(NDTreeNode *node, dimension_t dim) +{ + assert(!ndtree_node_is_leaf(node)); + // Compute the bounds of the children. + for (uint8_t i = 0; i < node->n; i++) { + NDTreeNode * child = node->children + i; + ndtree_node_compute_bounds_of_leaf(child, dim); + } +} + +/** + For a non-leaf node, whose children are all leaves with correct bounds, + (re-)compute the bounds of the node. This function does not + propagate the bounds up the tree. See ndtree_node_propagate_bounds_up() +*/ +static void +ndtree_node_compute_bounds_from_children(NDTreeNode *node, dimension_t dim) +{ + assert(!ndtree_node_is_leaf(node)); + // Initialize with the bounds of the first child. + NDTreeNode * child = node->children; + memcpy(node->ideal, child->ideal, dim * sizeof(double)); + memcpy(node->nadir, child->nadir, dim * sizeof(double)); + memcpy(node->midpoint, child->midpoint, dim * sizeof(double)); + + // Examine the bounds of the other children. + bool any_changed = false; + for (uint8_t i = 1; i < node->n; i++) { + child = node->children + i; + any_changed |= update_ideal(node, child->ideal, dim); + any_changed |= update_nadir(node, child->nadir, dim); + } + if (any_changed) + update_midpoint(node, dim); +} + + +/** + Copy bucket to the next element in node->bucket. It does NOT update bounds. +*/ +static inline void +bucket_append_move(NDTreeNode * restrict node, const SolutionsList * restrict bucket) +{ + assert(node->n < node->cap); + solutions_list_copy(&node->bucket[node->n], bucket); + node->n++; +} + +/** + bucket[] has max_bucket_size + 1 capacity, so we should never go over this + capacity and there is no need to reallocate. + + return true if it succeeds, false if not enough memory. + */ +static inline bool +bucket_append(NDTreeNode * restrict node, const double * restrict z, const void * restrict x, dimension_t dim) +{ + assert(node->n < node->cap); + if (!solutions_list_init(&node->bucket[node->n], z, x, dim)) + return false; + node->n++; + return true; +} + +static size_t ndtree_assert_valid_down(const NDTreeNode * node, dimension_t dim); +static void ndtree_assert_valid_up(const NDTreeNode *node, dimension_t dim); +static void ndtree_assert_valid_tree(const NDTreeArchive * tree); + +/** + Initialize a new leaf node. +*/ +static NDTreeNode * +ndtree_node_init(NDTreeNode *node, dimension_t dim, uint8_t max_bucket_size, NDTreeNode *parent) +{ + node->parent = parent; + // Allocate a single block for ideal/nadir/midpoint. + node->ideal = malloc(3U * dim * sizeof *node->ideal); + if (node->ideal == NULL) { + // PyErr_NoMemory(); + return NULL; + } + node->nadir = node->ideal + dim; + node->midpoint = node->nadir + dim; + node->coverage = 0; + // cap is max_bucket_size + 1 to be able to detect that the bucket is full. + if (ndtree_bucket_new(node, max_bucket_size + 1) == NULL) { + free(node->ideal); + return NULL; + } + return node; +} + +/** + Allocate and initialize a new leaf node. +*/ +static inline NDTreeNode * +ndtree_node_new(dimension_t dim, uint8_t max_bucket_size, NDTreeNode *parent) +{ + NDTreeNode *node = malloc(sizeof *node); + if (node == NULL) { + // PyErr_NoMemory(); + return NULL; + } + return ndtree_node_init(node, dim, max_bucket_size, parent); +} + +static inline NDTreeNode * +ndtree_init(NDTreeArchive *tree, dimension_t dim, + uint8_t max_children, uint8_t max_bucket_size, bool allow_duplicates) +{ + if (dim < 2 || dim >= MOOCORE_HV_DIMENSION_MAX) { + // PyErr_SetString(PyExc_ValueError, "number_of_objectives must be at least 2 and less than 255"); + return NULL; + } + + if (max_children == 0) + max_children = dim + 1; // default value + + if (max_bucket_size < max_children) + max_bucket_size = (uint8_t) MAX((int) max_children, 20); // default value + + tree->dim = dim; + tree->config.max_children = max_children; + tree->config.max_bucket_size = max_bucket_size; + tree->config.allow_duplicates = allow_duplicates; + tree->x_is_null = -1; // Don't know yet. + tree->root = ndtree_node_new(dim, max_bucket_size, NULL); + tree->revision = 0; + return tree->root; +} + +static inline NDTreeArchive * +ndtree_new(dimension_t dim, uint8_t max_children, uint8_t max_bucket_size, + bool allow_duplicates) +{ + NDTreeArchive * tree = malloc(sizeof(*tree)); + if (tree == NULL) { + free(tree); + return NULL; + } + if (!ndtree_init(tree, dim, max_children, max_bucket_size, allow_duplicates)) { + free(tree); + return NULL; + } + return tree; +} + +static inline bool +ndtree_init_with_single_solution(NDTreeArchive *restrict tree, + const double * restrict z, const void * restrict x, + dimension_t dim, const NDTreeConfig * restrict config) +{ + ndtree_init(tree, dim, config->max_children, config->max_bucket_size, config->allow_duplicates); + tree->x_is_null = (x == NULL) ? 1 : 0; + if (!bucket_append(tree->root, z, x, dim)) + return false; + assert(tree->root->n == 1); + ndtree_node_set_bounds(tree->root, z, dim); + DEBUG1(ndtree_assert_valid_tree(tree)); + return true; +} + +static void +ndtree_node_free(NDTreeNode *node); + +static void +ndtree_node_free_children(NDTreeNode *node) +{ + assert(node->children != NULL); + for (size_t i = 0; i < node->n; i++) + ndtree_node_free(&node->children[i]); + free(node->children); + node->children = NULL; +} +/** + Free the memory consumed by node and all its children. + +*/ +static void +ndtree_node_free(NDTreeNode *node) +{ + if (node == NULL) + return; + + free(node->ideal); + node->ideal = NULL; + node->cap = 0; + if (ndtree_node_is_leaf(node)) { + ndtree_bucket_free(node); + return; + } + ndtree_node_free_children(node); +} + +static void +ndtree_free_nodes(NDTreeArchive *tree) +{ + ndtree_node_free(tree->root); + free(tree->root); +} + +static inline bool +ndtree_node_is_empty(const NDTreeNode *node) +{ + assert(node != NULL); + return node->n == 0; +} + +static size_t +ndtree_node_coverage(const NDTreeNode *node) +{ + if (!ndtree_node_is_leaf(node)) + return node->coverage; + + size_t total = 0; + for (size_t i = 0; i < node->n; i++) + total += node->bucket[i].x.n; + return total; +} + +static size_t +ndtree_total_size(const NDTreeArchive *tree) +{ + return ndtree_node_coverage(tree->root); +} + +// Should we store in the tree the unique counts as well? +static size_t +ndtree_node_unique_count(const NDTreeNode *node) +{ + if (!ndtree_node_is_leaf(node)) { + size_t total = 0; + for (size_t i = 0; i < node->n; i++) + total += ndtree_node_unique_count(&node->children[i]); + return total; + } + return node->n; +} + +static size_t +ndtree_unique_size(const NDTreeArchive *tree) +{ + if (!tree->config.allow_duplicates) + return ndtree_total_size(tree); + return ndtree_node_unique_count(tree->root); +} + + +static void +ndtree_node_count_leaves(const NDTreeNode * restrict node, + size_t * restrict total, size_t * restrict count) +{ + if (ndtree_node_is_leaf(node)) { + *total += node->n; + (*count)++; + return; + } + + for (size_t i = 0; i < node->n; i++) + ndtree_node_count_leaves(node->children + i, total, count); +} + +static void +ndtree_node_count_internal(const NDTreeNode * restrict node, + size_t * restrict total, size_t * restrict count) +{ + if (!ndtree_node_is_leaf(node)) { + *total += node->n; + (*count)++; + for (size_t i = 0; i < node->n; i++) + ndtree_node_count_internal(node->children + i, total, count); + } +} + +static size_t +ndtree_node_max_height(const NDTreeNode *node) +{ + if (!ndtree_node_is_leaf(node)) { + size_t max_h = 0; + assert(node->n > 0); + for (size_t i = 0; i < node->n; i++) + max_h = MAX(max_h, ndtree_node_max_height(&node->children[i])); + return max_h + 1; + } + return 1; +} + +static size_t +ndtree_node_min_height(const NDTreeNode *node) +{ + if (!ndtree_node_is_leaf(node)) { + size_t min_h = SIZE_MAX; + assert(node->n > 0); + for (size_t i = 0; i < node->n; i++) + min_h = MIN(min_h, ndtree_node_min_height(&node->children[i])); + return min_h + 1; + } + return 1; +} + +static inline void +ndtree_print_stats(FILE * stream, const NDTreeArchive *tree) +{ + size_t total_size = ndtree_total_size(tree); + size_t unique_size = ndtree_unique_size(tree); + fprintf(stream, "# Min height = %zu\n", ndtree_node_min_height(tree->root)); + fprintf(stream, "# Max height = %zu\n", ndtree_node_max_height(tree->root)); + fprintf(stream, "# Total size = %zu\n", total_size); + fprintf(stream, "# Unique size = %zu\n", unique_size); + size_t total = 0, count = 0; + ndtree_node_count_leaves(tree->root, &total, &count); + assert(total == unique_size); + fprintf(stream, "# Number of leaves = %zu\n", count); + fprintf(stream, "# Occupancy = %g%%\n", + count > 0 ? (100. * (double)total / (double)(count * tree->config.max_bucket_size)) : 0); + total = 0; + count = 0; + ndtree_node_count_internal(tree->root, &total, &count); + fprintf(stream, "# Total children = %zu\n", total); + fprintf(stream, "# Number of internal nodes = %zu\n", count); + fprintf(stream, "# Occupancy = %g%%\n", + count > 0 ? (100. * (double)total / (double)(count * tree->config.max_children)) : 0); +} + +static inline void +ndtree_free(NDTreeArchive *tree) +{ + if (tree == NULL) + return; + ndtree_free_nodes(tree); + free(tree); +} + +static inline bool +ndtree_has_x_values(const NDTreeArchive *t) +{ + // -1 or 1 means that there are not x values to retrieve. + return t->x_is_null == 0; +} + +static void +ndtree_go_up_reducing_size(NDTreeNode *node, uint32_t removed) +{ + // Iterative instead of recursive because it is faster and helps inlining. + node = node->parent; + while (node) { + ASSUME(node->coverage >= removed); + node->coverage -= removed; + node = node->parent; + } +} + +static void +ndtree_go_up_increasing_size(NDTreeNode *node, uint32_t added) +{ + // Iterative instead of recursive because it is faster and helps inlining. + node = node->parent; + while (node) { + node->coverage += added; + node = node->parent; + } +} + + +static bool +ndtree_z_is_inside_box(const NDTreeNode *node, const double * restrict z, dimension_t dim) +{ + ASSUME(dim >= 2); + // FIXME: splitting the loop may help vectorization. + for (dimension_t i = 0; i < dim; i++) { + if (z[i] < node->ideal[i] || z[i] > node->nadir[i]) + return false; + } + return true; +} + +/** + If z is already in one of the buckets, append to the corresponding solution + list. Otherwise, reject it. + */ +static inline int +ndtree_node_bucket_append_if_present(NDTreeNode * restrict node, + const double * restrict z, const void * restrict x, + dimension_t dim) +{ + assert(ndtree_node_is_leaf(node)); + for (int i = 0; i < node->n; i++) { + if (all_equal_double(node->bucket[i].z, z, dim)) { + if (!solutions_list_append(&node->bucket[i], x)) + return ARCHIVE_INSERT_MEMORY_ERROR; + return ARCHIVE_INSERT_DUPLICATED; + } + } + return ARCHIVE_INSERT_REJECTED; +} + +/** + Append a new bucket (z, x) to node. The caller must test wrong_x_is_null(). + + This function does NOT update the bounds. The caller must handle that. +*/ +static int +ndtree_node_bucket_append(NDTreeNode * restrict node, + const double * restrict z, const void * restrict x, + dimension_t dim) +{ + assert(ndtree_node_is_leaf(node)); + assert(ndtree_node_bucket_append_if_present(node, z, x, dim) == ARCHIVE_INSERT_REJECTED); + if (!bucket_append(node, z, x, dim)) + return ARCHIVE_INSERT_MEMORY_ERROR; + return ARCHIVE_INSERT_ACCEPTED; +} + +static const NDTreeNode * +ndtree_find_node_of_vector(const NDTreeNode * restrict node, + const double * restrict z, dimension_t dim, + uint8_t * restrict bucket_index) +{ + ASSUME(dim >= 2); + if (ndtree_node_is_leaf(node)) { + for (uint8_t i = 0; i < node->n; i++) { + if (all_equal_double(node->bucket[i].z, z, dim)) { + *bucket_index = i; + return node; + } + } + return NULL; + } + for (size_t i = 0; i < node->n; i++) { + if (ndtree_z_is_inside_box(&node->children[i], z, dim)) { + const NDTreeNode * res = ndtree_find_node_of_vector(&node->children[i], z, dim, bucket_index); + if (res) + return res; + } + } + return NULL; +} + +static inline VoidList * +ndtree_find_solutions_of_vector(const NDTreeArchive * restrict tree, + const double * restrict z) +{ + uint8_t bucket_index; + const NDTreeNode *node = ndtree_find_node_of_vector( + tree->root, z, tree->dim, &bucket_index); + if (!node) + return NULL; + return &node->bucket[bucket_index].x; +} + +static inline bool +ndtree_find_exact_vector(const NDTreeArchive * restrict tree, + const double * restrict z) +{ + uint8_t bucket_index; + return ndtree_find_node_of_vector(tree->root, z, tree->dim, &bucket_index) != NULL; +} + +static bool +ndtree_node_dominates(const NDTreeNode * restrict node, + const double * restrict z, dimension_t dim) +{ + if (dominates(node->nadir, z, dim)) + return true; + if (weakly_dominates(z, node->ideal, dim)) + return false; + + if (dominates(node->ideal, z, dim) || dominates(z, node->nadir, dim)) { + if (ndtree_node_is_leaf(node)) { + for (size_t i = 0; i < node->n; i++) { + const double *bucket_z = node->bucket[i].z; + if (dominates(bucket_z, z, dim)) + return true; + if (weakly_dominates(z, bucket_z, dim)) + return false; + } + } else { + for (size_t i = 0; i < node->n; i++) { + if (ndtree_node_dominates(&node->children[i], z, dim)) + return true; + } + } + } + return false; +} + +static inline bool +ndtree_dominates(const NDTreeArchive * restrict tree, const double * restrict z) +{ + return !ndtree_node_is_empty(tree->root) && ndtree_node_dominates(tree->root, z, tree->dim); +} + +/** + Returns true if z dominates every point in the tree, false otherwise or if z + exists in the tree or the tree is empty. + +*/ +static inline bool +ndtree_dominated_by(const NDTreeArchive * restrict tree, const double * restrict z) +{ + if (ndtree_node_is_empty(tree->root)) + return false; + + const dimension_t dim = tree->dim; + const double * ideal = tree->root->ideal; + + /* The ideal at the root is tight as long as we never erase a nondominated + solution. Further down the tree a node's ideal may need tightening, as + an element leading to the ideal bound value may be dominated and + removed, but the dominator put in a different subtree, so the node with + the deletion could have its ideal now over-estimated. The root does not + have such a risk. */ + if (!weakly_dominates(z, ideal, dim)) + return false; + + return !tree->config.allow_duplicates + || !all_equal_double(z, ideal, dim) + || !ndtree_find_exact_vector(tree, z); +} + +static void +ndtree_remove_child_at(NDTreeNode * restrict parent, uint8_t index) +{ + NDTreeNode * removed _attr_maybe_unused = &parent->children[index]; + assert(ndtree_node_is_leaf(removed)); + // The child must be already empty. + assert(removed->bucket == NULL && removed->children == NULL && removed->ideal == NULL); + + parent->n--; +#if 0 + for (uint8_t i = index; i < parent->n; i++) { + parent->children[i] = parent->children[i+1]; + NDTreeNode * child = &parent->children[i]; + // We have to relink the grand children to their new parent. + if (!ndtree_node_is_leaf(child)) { + for (uint8_t k = 0; k < child->n; k++) { + child->children[k].parent = child; + } + } + } +#else + if (index < parent->n) { + // The order of the children is irrelevant, so overwrite with the last one. + parent->children[index] = parent->children[parent->n]; + NDTreeNode * child = &parent->children[index]; + // We have to relink the grand children to their new parent. + if (!ndtree_node_is_leaf(child)) { + for (uint8_t k = 0; k < child->n; k++) { + child->children[k].parent = child; + } + } + } +#endif +} + +/** + Changes the node type after prunning some children. +*/ +static NDTreeNode * +ndtree_remove_node_after_pruning(NDTreeNode *node, uint8_t max_bucket_size) +{ + if (node->n == 1) { // Replace tree with the only child. + NDTreeNode *child = node->children; + free(node->ideal); + child->parent = node->parent; + assert(node->bucket == NULL); + memcpy(node, child, sizeof(*node)); + free(child); + // Update the children's parents. + if (!ndtree_node_is_leaf(node)) { + for (uint8_t i = 0; i < node->n; i++) + node->children[i].parent = node; + } + assert(node->bucket == NULL || node->children == NULL); + } else if (node->n == 0) { // Becomes an empty leaf node. + node->coverage = 0; + ndtree_node_free_children(node); + assert(node->bucket == NULL); + assert(node->ideal != NULL); + if (node->parent != NULL) { + free(node->ideal); + node->ideal = NULL; + } + // If this is the root, we allocate an empty bucket. Otherwise, this + // node will be removed by ndtree_remove_child_at() later. + if (node->parent == NULL + && !ndtree_bucket_new(node, max_bucket_size + 1)) + return NULL; + } + return node; +} + +static NDTreeNode * +ndtree_closest_child(NDTreeNode *children, uint8_t n, const double * z, dimension_t dim) +{ + ASSUME(n < 255); + uint8_t which_inside[255]; + uint8_t num_inside = 0; + for (uint8_t i = 0; i < n; i++) { + NDTreeNode *child = &children[i]; + assert(child->n > 0); + if (ndtree_z_is_inside_box(child, z, dim)) { + which_inside[num_inside] = i; + num_inside++; + } + } + + if (num_inside > 0) { + uint8_t i = 0; + NDTreeNode * closest_child = &children[which_inside[i]]; + if (num_inside > 1) { + double closest_distance = squared_distance(closest_child->midpoint, z, dim); + for (i++; i < num_inside; i++) { + NDTreeNode *child = &children[which_inside[i]]; + double distance = squared_distance(child->midpoint, z, dim); + if (distance < closest_distance) { + closest_distance = distance; + closest_child = child; + } + } + } + return closest_child; + } + + // z is not inside any box, so we have to calculate all distances. + NDTreeNode * closest_child = &children[0]; + double closest_distance = squared_distance(closest_child->midpoint, z, dim); + for (uint8_t i = 1; i < n; i++) { + NDTreeNode *child = &children[i]; + double distance = squared_distance(child->midpoint, z, dim); + if (distance < closest_distance) { + closest_distance = distance; + closest_child = child; + } + } + return closest_child; +} + + + +/** + Splits a leaf node by converting it to an internal node and distributing the + bucket elements among newly created children. + + Complexity O(max(n_children * n_children, max_bucket_size) * max_bucket_size + ?) + */ +static int +ndtree_node_split(NDTreeNode *node, dimension_t dim, + uint8_t max_bucket_size, uint8_t max_children) +{ + ASSUME(node->n < 255); + /* At this point, the bounds may be wrong because we inserted a new z and + it triggered the split, but we did not update the bounds. */ + uint8_t bucket_n = node->n; + ASSUME(bucket_n >= max_children); + assert(ndtree_node_is_leaf(node)); + ASSUME(node->coverage == 0); + for (uint8_t i = 0; i < bucket_n; i++) + node->coverage += node->bucket[i].x.n; + ASSUME(node->coverage > 0); + + // distance[i]: sum of distances of solution i to all other solutions. + double *distance_sum = calloc(bucket_n + bucket_n * bucket_n, sizeof(*distance_sum)); + double *distance_matrix = distance_sum + bucket_n; + NDTreeNode * children = calloc(max_children, sizeof *children); + if (distance_sum == NULL || children == NULL) { + free(distance_sum); + free(children); + return ARCHIVE_INSERT_MEMORY_ERROR; + + } + const SolutionsList *bucket = node->bucket; + /* FIXME: We could probably vectorize this loop by copying the midpoints to + a temporary buffer, then computing all distances at once. But is the + copy worth the potential speed-up? */ + for (unsigned i = 0; i < bucket_n; i++) { + for (unsigned j = i + 1; j < bucket_n; j++) { + double d = squared_distance(bucket[i].z, bucket[j].z, dim); + distance_matrix[i * bucket_n + j] = d; + distance_matrix[j * bucket_n + i] = d; + distance_sum[i] += d; + distance_sum[j] += d; + } + } + + NDTreeNode * child = &children[0]; + if (!ndtree_node_init(child, dim, max_bucket_size, node)) { + free(distance_sum); + free(children); + return ARCHIVE_INSERT_MEMORY_ERROR; + } + uint8_t seed_rep = (uint8_t) which_max(distance_sum, bucket_n); + ASSUME(seed_rep < bucket_n); + // We do not update bounds until the end of this function, when we know the + // child assigned to each bucket. + bucket_append_move(child, &bucket[seed_rep]); + + // Store the buckets assigned to children and the ones not yet assigned. + uint8_t bucket_of_child[255]; + // We store all buckets indexes, then swap to the front the ones assigned + // to a child. + for (uint8_t i = 0; i < bucket_n; i++) + bucket_of_child[i] = i; + bucket_of_child[0] = seed_rep; + bucket_of_child[seed_rep] = 0; + + /* FIXME: It may be faster to find the two buckets closest to each other, + then put them together, then treat them as a single bucket, and keep + doing that until n_buckets == max_children, then simply assign each + (virtual) bucket to each child. Then update bounds in one go. */ + for (unsigned i = 1; i < max_children; i++) { + double max_distance = -1; + size_t seed_rep = SIZE_MAX; + for (unsigned k = i; k < bucket_n; k++) { + uint8_t bucket_k = bucket_of_child[k]; + double accumulator = 0.0; + for (unsigned j = 0; j < i; j++) + accumulator += distance_matrix[bucket_k * bucket_n + bucket_of_child[j]]; + + if (accumulator > max_distance) { + max_distance = accumulator; + seed_rep = k; + } + } + ASSUME(seed_rep < SIZE_MAX); // There must be one not added yet. + child = &children[i]; + if (!ndtree_node_init(child, dim, max_bucket_size, node)) { + free(distance_sum); + free(children); + return ARCHIVE_INSERT_MEMORY_ERROR; + } + + bucket_append_move(child, &bucket[bucket_of_child[seed_rep]]); + if (i != seed_rep) + SWAP(bucket_of_child[i], bucket_of_child[seed_rep]); + } + /* At this point all children have one bucket, given by + bucket_of_child[i] when i < max_children, so we can re-use the + previously calculated distances. */ + if (max_children < bucket_n) { + uint8_t bucket_k = bucket_of_child[max_children]; + uint8_t min_child = 0; + double min_dist = distance_matrix[min_child * bucket_n + bucket_k]; + for (uint8_t j = 1; j < max_children; j++) { + double dist_j = distance_matrix[bucket_of_child[j] * bucket_n + bucket_k]; + if (dist_j < min_dist) { + min_child = j; + min_dist = dist_j; + } + } + child = &children[min_child]; + bucket_append_move(child, &bucket[bucket_k]); + // We do not need to SWAP. + } + + free(distance_sum); + node->children = children; + node->n = max_children; + node->cap = max_children; + ndtree_node_compute_bounds_of_children(node, dim); + if (max_children < bucket_n - 1) { + /* We are not done yet! The are still buckets not assigned to children, + but now we need to update bounds every time we add a new bucket. */ + + /* Now some children have two buckets, so we have to rely on + ndtree_closest_child(), which is quite wasteful for children that + have only one. */ + for (uint8_t k = max_children + 1; k < bucket_n; k++) { + uint8_t bucket_k = bucket_of_child[k]; + const double *z = bucket[bucket_k].z; + NDTreeNode *target = ndtree_closest_child(children, max_children, z, dim); + bucket_append_move(target, &bucket[bucket_k]); + bool ideal_changed = true, nadir_changed = true; + ndtree_node_update_bounds(target, z, dim, &ideal_changed, &nadir_changed); + } + } + ndtree_node_compute_bounds_from_children(node, dim); + ndtree_node_propagate_bounds_up(node, dim); + // FIXME: Ideally, we would reuse the memory. + free(node->bucket); + node->bucket = NULL; + DEBUG1(ndtree_assert_valid_up(node, dim)); + DEBUG1(ndtree_assert_valid_down(node, dim)); + return ARCHIVE_INSERT_ACCEPTED; +} + +/** + Insert entry in node. If node is not a leaf node, try again with the closest + child and so on. Once it is inserted in a leaf node, maybe split it if it + has grown too large. + +*/ +static int +ndtree_node_insert(NDTreeArchive * restrict tree, NDTreeNode * restrict node, + const double * restrict z, const void * restrict x) +{ + const dimension_t dim = tree->dim; + const uint8_t max_bucket_size = tree->config.max_bucket_size; + while (!ndtree_node_is_leaf(node)) { + node->coverage++; + node = ndtree_closest_child(node->children, node->n, z, dim); + } + ASSUME(tree->x_is_null == 1 || tree->x_is_null == 0); + int res = wrong_x_is_null(&tree->x_is_null, x); + if (res < 0) + return res; + + res = ndtree_node_bucket_append(node, z, x, dim); + assert(res == ARCHIVE_INSERT_ACCEPTED); + ASSUME(ndtree_node_is_leaf(node)); + // Only update the bounds if we do not split the node. + if (node->n == 1) { + // The node was empty before z + // This cannot happen in the current code, because we never add to an empty child. + assert(false); + ndtree_node_set_bounds(node, z, dim); + ndtree_node_propagate_bounds_up(node, dim); + DEBUG2(ndtree_assert_valid_tree(tree)); + return res; + } else if (node->n <= max_bucket_size) { + ndtree_node_update_bounds_up(node, z, dim); + DEBUG2(ndtree_assert_valid_tree(tree)); + return res; + } + return ndtree_node_split(node, dim, max_bucket_size, tree->config.max_children); +} + +/** + Clear the subtree that has node as a parent and convert node to an empty + leaf node. If node is the root of the tree, the tree is reset to its empty + state. + + This does not free node nor its bucket, but it does free all solutions + lists. +*/ +static bool +ndtree_node_clear_to_empty_leaf(NDTreeNode *node, uint8_t max_bucket_size) +{ + // FIXME: This function should record what is being removed, so we can + // return it. + if (ndtree_node_is_leaf(node)) { + /* An empty leaf has ->n=0, but ->bucket is already initialized to + ->cap size and each ->bucket[i].x.n = 0 with ->bucket[i].x.list == NULL. */ + assert(node->bucket != NULL); + /* FIXME: The solutions lists may be empty if we moved these buckets + elsewhere. How can we detect that? */ + for (int i = 0; i < node->n; i++) + solutions_list_free(&node->bucket[i]); + + } else { + ndtree_node_free_children(node); + if (!ndtree_bucket_new(node, max_bucket_size + 1)) + return false; + } + node->n = 0; + node->coverage = 0; + return true; +} + +static uint32_t +ndtree_append_subtree_to_list(FlexBucket ** restrict list, NDTreeNode * restrict node) +{ + if (ndtree_node_is_leaf(node)) + return (uint32_t) flex_bucket_append_move(list, node->bucket, node->n); + + uint32_t total = 0; + for (unsigned i = 0; i < node->n && list; i++) + total += ndtree_append_subtree_to_list(list, &node->children[i]); + return total; +} + +static int +ndtree_add_duplicate(signed char x_is_null, NDTreeNode * restrict node, + uint8_t bucket_index, const void * restrict x) +{ + ASSUME(x_is_null == 1 || x_is_null == 0); + int res = wrong_x_is_null(&x_is_null, x); + if (res < 0) + return res; + if (!solutions_list_append(&node->bucket[bucket_index], x)) + return ARCHIVE_INSERT_MEMORY_ERROR; + ndtree_go_up_increasing_size(node, 1); + return ARCHIVE_INSERT_DUPLICATED; +} + +static int +ndtree_update_node(NDTreeArchive * restrict tree, NDTreeNode * restrict node, + const double * restrict z, const void * restrict x, + bool check, bool allow_duplicates, FlexBucket **displaced) +{ + const dimension_t dim = tree->dim; + if (check && weakly_dominates(node->nadir, z, dim)) { + // Maybe there is one solution and it is a duplicate. + if (allow_duplicates && ndtree_node_is_leaf(node) && node->n == 1 + && all_equal_double(z, node->bucket[0].z, dim)) { + return ndtree_add_duplicate(tree->x_is_null, node, 0, x); + } + return ARCHIVE_INSERT_REJECTED; + } + + bool clears_subtree = allow_duplicates + ? dominates(z, node->ideal, dim) + : weakly_dominates(z, node->ideal, dim); + if (clears_subtree) { + uint32_t removed = ndtree_append_subtree_to_list(displaced, node); + if (!*displaced) + return ARCHIVE_INSERT_MEMORY_ERROR; + + if (node->parent != NULL) { + ndtree_go_up_reducing_size(node, removed); + // Free bucket/children and subtrees. Mark the node as + // uninitialized. The caller will call ndtree_remove_child_at() to + // remove it from node->parent->children[] + ndtree_node_free(node); + node->coverage = 0; + node->n = 0; + DEBUG2(ndtree_assert_valid_tree(tree)); + } else { + // This completely clears the tree up to an empty root leaf. + /* FIXME: Why not reuse the subtree for inserting the new entry OR + detach the node from the tree and add it to displaced? */ + if (!ndtree_node_clear_to_empty_leaf(node, tree->config.max_bucket_size)) { + return ARCHIVE_INSERT_MEMORY_ERROR; + } + DEBUG2(ndtree_assert_valid_tree(tree)); + } + return ARCHIVE_INSERT_ACCEPTED; + } + + DEBUG2(ndtree_assert_valid_tree(tree)); + bool z_wdom_nadir = weakly_dominates(z, node->nadir, dim); + bool ideal_wdom_z = weakly_dominates(node->ideal, z, dim); + if (!ideal_wdom_z && !z_wdom_nadir) + return ARCHIVE_INSERT_ACCEPTED; + + // If !ideal_wdom_z, then nothing in this node can weakly dominate z, so do + // not check. + check = check && ideal_wdom_z; + if (ndtree_node_is_leaf(node)) { + if (z_wdom_nadir) { + uint32_t removed = 0; + for (uint8_t i = 0; i < node->n; ) { + SolutionsList *bucket = &node->bucket[i]; + if (weakly_dominates(z, bucket->z, dim)) { + if (all_equal_double(z, bucket->z, dim)) { + if (allow_duplicates) + return ndtree_add_duplicate(tree->x_is_null, node, i, x); + // equal not allowed. + return ARCHIVE_INSERT_REJECTED; + } + // z dominates bucket_z + check = false; // No other bucket can dominate z. + removed += bucket->x.n; + flex_bucket_append_move(displaced, bucket, 1); + if (!*displaced) + return ARCHIVE_INSERT_MEMORY_ERROR; + node->n--; + if (i != node->n) + SWAP(node->bucket[i], node->bucket[node->n]); + // Do not increment i. + } else if (check && weakly_dominates(bucket->z, z, dim)) { + return ARCHIVE_INSERT_REJECTED; + } else { + i++; + } + } + if (removed > 0) { + ndtree_go_up_reducing_size(node, removed); + if (node->n == 0 && node->parent != NULL) { + /* Since node->n == 0, the ndtree_node_free() will not call + solutions_list_free(), which is what we want because they + are now in *displaced. */ + ndtree_node_free(node); + } + } + } else if (check) { + // If !z_wdom_nadir then z cannot dominated anything in this leaf, + // but it could be dominated. + for (uint8_t i = 0; i < node->n; i++) { + if (weakly_dominates(node->bucket[i].z, z, dim)) + return ARCHIVE_INSERT_REJECTED; + } + } + } else { // NOT a leaf node. + uint8_t index = 0; + while (index < node->n) { + int res = ndtree_update_node(tree, &node->children[index], z, x, + check, allow_duplicates, displaced); + /* We may have a slightly invalid tree here because we may have + removed all children, so coverage = 0 but n = 1. That is, this + looks like a leaf but it is just empty. */ + if (res != ARCHIVE_INSERT_ACCEPTED) + return res; + + if (ndtree_node_is_empty(&node->children[index])) { + ndtree_remove_child_at(node, index); + DEBUG2(ndtree_assert_valid_tree(tree)); + } else { + index++; + assert(index < 255); + DEBUG2(ndtree_assert_valid_tree(tree)); + } + } + if (ndtree_remove_node_after_pruning(node, tree->config.max_bucket_size) == NULL) + return ARCHIVE_INSERT_MEMORY_ERROR; + /* At this point, we may have have converted node to a (possibly empty) + leaf node. */ + DEBUG2(ndtree_assert_valid_tree(tree)); + } + return ARCHIVE_INSERT_ACCEPTED; +} + +static int +ndtree_add_impl(NDTreeArchive * restrict tree, + const double * restrict z, const void * restrict x, + bool check, bool allow_duplicates, FlexBucket **displaced) +{ + const dimension_t dim = tree->dim; + DEBUG1(ndtree_assert_valid_tree(tree)); + + if (ndtree_node_is_empty(tree->root)) { + if (tree->x_is_null == -1) + tree->x_is_null = (x == NULL); + int res = ndtree_node_bucket_append(tree->root, z, x, dim); + assert(res < 0 || res == ARCHIVE_INSERT_ACCEPTED); + assert(tree->root->n == 1); + ndtree_node_set_bounds(tree->root, z, dim); + DEBUG1(ndtree_assert_valid_tree(tree)); + // At this point we have modified the tree. + tree->revision++; + return res; + } + + int res = ndtree_update_node(tree, tree->root, z, x, check, allow_duplicates, displaced); + if (res != ARCHIVE_INSERT_ACCEPTED) + return res; + // At this point we have modified the tree. + tree->revision++; + // It is empty only if we have removed everything. + if (ndtree_node_is_empty(tree->root)) { + res = wrong_x_is_null(&tree->x_is_null, x); + if (res < 0) + return res; + res = ndtree_node_bucket_append(tree->root, z, x, dim); + assert(tree->root->n == 1); + ndtree_node_set_bounds(tree->root, z, dim); + DEBUG1(ndtree_assert_valid_tree(tree)); + return res; + } + // FIXME: Now we have to traverse the tree again to see where to insert + // the node. Why not try to insert it in update_node() ? Or at least + // return which subtree is the most promising and start from that. + DEBUG1(ndtree_assert_valid_tree(tree)); + res = ndtree_node_insert(tree, tree->root, z, x); + DEBUG1(ndtree_assert_valid_tree(tree)); + return res; +} + +/** + Insert new solution (z,x) and get the displaced solutions. + + *displaced must be NULL or already properly initialized. + + Returns one of the values of archive_insert_result_t. +*/ +static inline int +ndtree_add(NDTreeArchive * restrict tree, + const double * restrict z, const void * restrict x, + bool check, FlexBucket **displaced) +{ + return ndtree_add_impl(tree, z, x, check, tree->config.allow_duplicates, displaced); +} + +static inline int +ndtree_add_discard(NDTreeArchive * restrict tree, + const double * restrict z, const void * restrict x) +{ + // FIXME: Avoid using displaced at all. + FlexBucket *displaced = flex_bucket_new(4); + int res = ndtree_add_impl(tree, z, x, /*check=*/true, tree->config.allow_duplicates, &displaced); + flex_bucket_free(displaced); + return res; +} + + +static int +ndtree_add_displaced_helper(NDTreeArchive * restrict tree, + const double * restrict z, FlexBucket **displaced) +{ + const void * x = NULL; + int res = ndtree_add_impl(tree, z, x, /*check=*/false, /*allow_duplicates=*/false, displaced); + assert(res == ARCHIVE_INSERT_ACCEPTED); + return res; +} + +/** + Insert incoming solutions into tree and return displaced nodes in subtree. + It does not allocate any new memory! +*/ +static inline int +ndtree_insert_displaced(NDTreeArchive *tree, FlexBucket **displaced_p, bool *all_displaced) +{ + assert(*displaced_p != NULL); + if (*displaced_p == NULL) + return ARCHIVE_INSERT_MEMORY_ERROR; + + FlexBucket * incoming = *displaced_p; + FlexBucket * displaced = NULL; + // FIXME: There must be a better way than adding the points one by one! + const bool tree_was_empty = ndtree_node_is_empty(tree->root); + _attr_maybe_unused size_t prev_tree_unique; + DEBUG1(prev_tree_unique = ndtree_unique_size(tree)); + + const bool incoming_x_is_null = (incoming->solutions[0].x.list == NULL); + // Failing this assert is a programming error and should never happen. + assert(tree->x_is_null == -1 || tree->x_is_null == incoming_x_is_null); + tree->x_is_null = 1; // Disable temporarily + for (size_t i = 0; i < incoming->n_unique; i++) { + // We should not have empty solutions. + assert(incoming->solutions[i].x.n > 0); + const double * z = incoming->solutions[i].z; + // FIXME: We should return the node where z is inserted and the bucket + // index, so we can avoid the ndtree_find_node_of_vector() call below. + // FIXME: We should also not allocate a copy of z in this case, so we + // do not need to free it below. + int res = ndtree_add_displaced_helper(tree, z, &displaced); + if (res < 0) + return res; + assert(res == ARCHIVE_INSERT_ACCEPTED); + + if (incoming_x_is_null && incoming->solutions[i].x.n == 1) { + free((void*)incoming->solutions[i].z); + continue; // We are done with this solutions[i]. + } + + /* Instead of inserting duplicates, find the just created SolutionsList + in the tree and modify it. */ + uint8_t bucket_index = -1; + // const_cast so we can update z_node. + NDTreeNode * z_node = (NDTreeNode *) ndtree_find_node_of_vector( + tree->root, z, tree->dim, &bucket_index); + SolutionsList * z_bucket = &z_node->bucket[bucket_index]; + free((void*)z_bucket->z); // allocated by ndtree_add_displaced_helper() + assert(z_bucket->x.list == NULL); // Did we allocate this list by mistake? + if (incoming_x_is_null) { + // We have duplicates but no x values, just copy z and the counter. + z_bucket->z = incoming->solutions[i].z; + z_bucket->x.n = incoming->solutions[i].x.n; + } else { + // We have duplicates and x values, just copy the whole thing. + solutions_list_copy(z_bucket, &incoming->solutions[i]); + } + // n - 1 because we already counted 1 when inserting z. + ndtree_go_up_increasing_size(z_node, z_bucket->x.n - 1); + } + tree->x_is_null = incoming_x_is_null; // Reset. + /* All displaced points will be inserted, so if the size of the tree is the + same as the number of displaced, all previous points in the tree have + been displaced. */ + assert(ndtree_unique_size(tree) >= incoming->n_unique); + assert(ndtree_unique_size(tree) > 0); + *all_displaced = !tree_was_empty && (ndtree_unique_size(tree) == incoming->n_unique); + assert(!*all_displaced || (displaced && prev_tree_unique == displaced->n_unique)); + assert(tree_was_empty || !displaced || displaced->n_unique <= prev_tree_unique); + + free(incoming); + *displaced_p = displaced; + DEBUG1(ndtree_assert_valid_tree(tree)); + return ARCHIVE_INSERT_ACCEPTED; +} + +static inline bool +ndtree_init_with_displaced(NDTreeArchive *tree, FlexBucket *displaced, + dimension_t dim, const NDTreeConfig * config) +{ + ndtree_init(tree, dim, config->max_children, config->max_bucket_size, config->allow_duplicates); + bool all_displaced; + int res = ndtree_insert_displaced(tree, &displaced, &all_displaced); + DEBUG1(ndtree_assert_valid_tree(tree)); + assert(res != ARCHIVE_INSERT_REJECTED && res != ARCHIVE_INSERT_DUPLICATED); + return res == ARCHIVE_INSERT_ACCEPTED; +} + +static inline size_t +ndtree_node_get_unique_vectors(const NDTreeNode *node, const double *** z_list, + dimension_t dim) +{ + if (ndtree_node_is_leaf(node)) { + for (size_t i = 0; i < node->n; i++) { + assert(node->bucket[i].z != NULL); + **z_list = node->bucket[i].z; + (*z_list)++; + } + return node->n; + } + size_t total = 0; + for (size_t i = 0; i < node->n; i++) + total += ndtree_node_get_unique_vectors(&node->children[i], z_list, dim); + return total; +} + +/** + z_list must be an array of pointers to 'const double' of length treap_unique_size(t). + + That is, z_list does not contain the actual data but just the pointers + stored in the tree. Since the tree does not stored the data either, the + data may get deleted by something else if not used or copied immediately. +*/ +static inline size_t +ndtree_get_unique_vectors(const NDTreeArchive *tree, const double ** z_list) +{ + return ndtree_node_get_unique_vectors(tree->root, &z_list, tree->dim); +} + + +static inline void +ndtree_node_get_x_values(const NDTreeNode *node, const void *** x_list) +{ + if (ndtree_node_is_leaf(node)) { // FIXME: Why check leafs first? + for (size_t i = 0; i < node->n; i++) + solutions_list_get_x_values(&node->bucket[i], x_list); + return; + } + + for (size_t i = 0; i < node->n; i++) + ndtree_node_get_x_values(&node->children[i], x_list); + +} + +static inline void +ndtree_get_x_values(const NDTreeArchive *arch, const void ** x_list) +{ + if (!ndtree_has_x_values(arch)) + return; + ndtree_node_get_x_values(arch->root, &x_list); +} + +static inline void +ndtree_node_get_contents(const NDTreeNode *node, + double ** z_list, const void *** x_list, + dimension_t dim) +{ + if (ndtree_node_is_leaf(node)) { // FIXME: Why check leafs first? + for (size_t i = 0; i < node->n; i++) + solutions_list_get_contents(&node->bucket[i], z_list, x_list, dim); + return; + } + + for (size_t i = 0; i < node->n; i++) + ndtree_node_get_contents(&node->children[i], z_list, x_list, dim); +} + +/** + z_list must be an array of double allocated with sufficient space. Vectors + are copied to this array. + + x_list must be of array of void *. Only the pointers are copied. +*/ +static inline void +ndtree_get_contents(const NDTreeArchive * restrict tree, double * z_list, + const void ** x_list) +{ + if (!ndtree_node_is_empty(tree->root)) { + ndtree_node_get_contents(tree->root, &z_list, + ndtree_has_x_values(tree) ? &x_list : NULL, + tree->dim); + } +} + + +static size_t +ndtree_assert_valid_down(const NDTreeNode * node, dimension_t dim) +{ + size_t total = 0; + if (ndtree_node_is_leaf(node)) { + for (size_t i = 0; i < node->n; i++) { + total += node->bucket[i].x.n; + for (dimension_t d = 0; d < dim; d++) { + assert(node->ideal[d] <= node->bucket[i].z[d]); + assert(node->nadir[d] >= node->bucket[i].z[d]); + } + } + return total; + } + for (size_t i = 0; i < node->n; i++) { + assert(node->children[i].parent == node); + total += ndtree_assert_valid_down(&node->children[i], dim); + for (dimension_t d = 0; d < dim; d++) { + assert(node->ideal[d] <= node->children[i].ideal[d]); + assert(node->nadir[d] >= node->children[i].nadir[d]); + } + } + assert(total == node->coverage); + return total; +} + +static void +ndtree_assert_valid_tree(const NDTreeArchive * tree) +{ + size_t total = ndtree_assert_valid_down(tree->root, tree->dim); + size_t unique_size = ndtree_unique_size(tree); + assert (total >= unique_size); + if (total == 0) + return; + + assert(tree->x_is_null != -1); + double * points = malloc(tree->dim * unique_size * sizeof(*points)); + const double ** z_list = malloc(unique_size * sizeof(*z_list)); + _attr_maybe_unused size_t new_unique = ndtree_get_unique_vectors(tree, z_list); + assert(unique_size == new_unique); + for (size_t i = 0; i < unique_size; i++) { + for (size_t d = 0; d < tree->dim; d++) { + assert(z_list[i] != NULL); + points[i * tree->dim + d] = z_list[i][d]; + } + } + const int * minmax = minmax_minimise(tree->dim); + _attr_maybe_unused size_t pos = + find_dominated_point_(points, unique_size, tree->dim, + /* keep_weakly=*/false, AGREE_NONE, minmax); + free((void*)minmax); + assert(pos >= unique_size); + free(points); + free(z_list); +} + +static void +ndtree_assert_valid_up(const NDTreeNode *node, dimension_t dim) +{ + while (node->parent != NULL) { + assert(!ndtree_node_is_leaf(node->parent)); + node = node->parent; + for (size_t i = 0; i < node->n; i++) { + assert(node->children[i].parent == node); + for (dimension_t d = 0; d < dim; d++) { + assert(node->ideal[d] <= node->children[i].ideal[d]); + assert(node->nadir[d] >= node->children[i].nadir[d]); + } + } + } +} + +#include "nd_tree_iterator.h" + +#endif // ND_TREE_ARCHIVE_H diff --git a/c/nd_tree_iterator.h b/c/nd_tree_iterator.h new file mode 100644 index 000000000..553121dd1 --- /dev/null +++ b/c/nd_tree_iterator.h @@ -0,0 +1,119 @@ +#ifndef NDTREE_ITERATOR_H +#define NDTREE_ITERATOR_H + +typedef struct NDTreeIterator { + const NDTreeNode *current_node; + size_t bucket_index; + size_t solution_index; + // To detect modifications of the tree while iterating. + uint64_t revision; + const NDTreeArchive *arch; +} NDTreeIterator; + +static void +ndtree_iter_init(NDTreeIterator * it, const NDTreeArchive *tree) +{ + it->bucket_index = 0; + it->solution_index = 0; + it->revision = tree->revision; + it->arch = tree; + // Find first leaf. + const NDTreeNode *node = tree->root; + if (node != NULL) { + while (!ndtree_node_is_leaf(node)) + node = &node->children[0]; + } + it->current_node = node; +} + +static inline NDTreeIterator * +ndtree_iter_new(const NDTreeArchive *tree) +{ + if (tree == NULL) + return NULL; + + NDTreeIterator *it = malloc(sizeof(*it)); + if (it == NULL) + return NULL; + + ndtree_iter_init(it, tree); + return it; +} + +static inline void +ndtree_iter_free(NDTreeIterator *it) +{ + free(it); +} + +static const NDTreeNode * +ndtree_next_leaf(const NDTreeNode *node) +{ + while (node != NULL) { + // Descend to the next subtree if there is a next sibling. + const NDTreeNode *parent = node->parent; + if (!parent) + return NULL; + + assert(node - parent->children >= 0); + assert(node - parent->children < 255); + uint8_t i = 1 + (uint8_t)(node - parent->children); + + assert(parent->n > 0); + if (i < parent->n) { + // Move to the next sibling, then descend to its leftmost leaf. + node = &parent->children[i]; + + while (!ndtree_node_is_leaf(node)) + node = &node->children[0]; + + return node; + } + // No more siblings here. Move up and look for one. + node = parent; + } + return NULL; +} + + +/** + 1 node returned + 0 end of iteration + -1 error +*/ +static inline int +ndtree_iter_next(NDTreeIterator *it, const double **z, const void **x) +{ + if (it == NULL || z == NULL || it->revision != it->arch->revision) + return -1; + + while (it->current_node != NULL) { + const NDTreeNode *node = it->current_node; + + assert(ndtree_node_is_leaf(node)); + // Emit solutions from the current leaf. + while (it->bucket_index < node->n) { + const SolutionsList *bucket = &node->bucket[it->bucket_index]; + if (it->solution_index < bucket->x.n) { + *z = bucket->z; + if (x) + *x = bucket->x.list[it->solution_index]; + it->solution_index++; + return 1; + } + it->bucket_index++; + it->solution_index = 0; + } + + // Current leaf exhausted. + it->current_node = ndtree_next_leaf(node); + it->bucket_index = 0; + it->solution_index = 0; + } + *z = NULL; + if (x) + *x = NULL; + return 0; +} + +#endif // NDTREE_ITERATOR_H diff --git a/c/nondominated.c b/c/nondominated.c index 112458daf..a947934fb 100644 --- a/c/nondominated.c +++ b/c/nondominated.c @@ -6,18 +6,7 @@ *****************************************************************************/ #include "nondominated.h" -typedef const double avl_item_t; -typedef struct avl_node_t { - struct avl_node_t *next; - struct avl_node_t *prev; - struct avl_node_t *parent; - struct avl_node_t *left; - struct avl_node_t *right; - avl_item_t *item; - unsigned char depth; -} avl_node_t; - -#include "avl_tiny.h" +#include "treap2d_fixed_size.h" /** 3D dimension-sweep algorithm by H. T. Kung, F. Luccio, and F. P. Preparata. @@ -38,73 +27,35 @@ find_nondominated_3d_impl_sorted(const double ** restrict rows, size_t size, const bool find_dominated) { ASSUME(size > 1); - /* FIXME: The AVL-tree is the bottleneck of this algorithm. A Treap - [R. Seidel and C. R. Aragon. Randomized search trees. Algorithmica, - 16:464–497, 1996] may be far more efficient by allowing to remove a - range of items without rebalancing. See - https://alexdremov.me/treap-algorithm-explained/ - Instead of using randomized priorities, hash the node pointer or the node index. - */ - avl_tree_t tree; - avl_init_tree(&tree, qsort_cmp_pdouble_asc_x_nonzero); - avl_node_t * tnodes = malloc((size+1) * sizeof(*tnodes)); - avl_node_t * node = tnodes; - node->item = rows[0]; - avl_insert_top(&tree, node); - - const double sentinel[] = { INFINITY, -INFINITY }; - (++node)->item = sentinel; - avl_insert_after(&tree, node - 1, node); + Treap2D * tree = treap2d_new(size); + assert(tree != NULL); + double pk0 = rows[0][0], pk1 = rows[0][1], pk2 = rows[0][2]; + treap2d_init_with_single_node(tree, pk0, pk1); // In this context, size means "no dominated solution found". size_t new_size = size; size_t prev_dominated = false; - double pk0 = rows[0][0], pk1 = rows[0][1], pk2 = rows[0][2]; for (size_t j = 1; j < size; j++) { const double * restrict pj = rows[j]; DEBUG2(printf_point("pj = [ ", pj, 3, " ], ")); const double pj0 = pj[0], pj1 = pj[1], pj2 = pj[2]; if ((pk0 > pj0) | (pk1 > pj1)) { // Check if pj is dominated by a point in the tree. - avl_node_t * nodeaux; - int res = avl_search_closest(&tree, pj, &nodeaux); - assert(res != 0); - if (res > 0 || nodeaux->prev) { - const double * restrict prev; - if (res > 0) { // nodeaux goes before pj - prev = nodeaux->item; - nodeaux = nodeaux->next; - DEBUG2(printf_point("res > 0: prev: ", prev, 3, "\n")); - } else { // nodeaux goes after pj, so move to the previous one. - prev = nodeaux->prev->item; - DEBUG2(printf_point("res < 0: prev: ", prev, 3, "\n")); - } - assert(prev[0] != sentinel[0]); - assert(prev[0] <= pj0); - if (prev[1] <= pj1) - goto j_is_dominated; - } + /* In a valid 2-D frontier, x increases and y decreases. Therefore + the only existing point that can dominate pj is the frontier + predecessor with the largest key <= pj0. */ + TreapNode *pred = treap2d_find_le(tree, pj0); + if (pred != NULL && treap2d_node_get_y(pred) <= pj1) + goto j_is_dominated; - // pj is NOT dominated by a point in the tree. - const double * restrict point = nodeaux->item; - assert(pj0 <= point[0]); - // Delete everything in the tree that is dominated by pj. - while (pj1 <= point[1]) { - DEBUG2(printf_point("delete point: ", point, 3, "\n")); - assert(pj0 <= point[0]); - nodeaux = nodeaux->next; - point = nodeaux->item; - /* FIXME: A possible speed up is to delete without rebalancing - the tree because avl_insert_before() will rebalance, but we - need to know which is the highest node that needs - rebalancing. */ - avl_unlink_node(&tree, nodeaux->prev); - } - DEBUG2((point == sentinel) - ? printf_point("insert before sentinel: ", sentinel, 2, "\n") - : printf_point("insert before point: ", point, 3, "\n")); - (++node)->item = pj; - avl_insert_before(&tree, nodeaux, node); + /* pj is not dominated by an existing frontier point. + + Insert it and detach every existing point dominated by it. + + The returned treap contains exactly those displaced nodes, but + we do not need to traverse it here because those nodes will + never again participate in dominance queries. */ + treap2d_insert(tree, pj0, pj1); // Fall-through to j_is_NOT_dominated. } // Handle duplicates and points that are dominated by the immediate previous one. else if (!keep_weakly // Don't keep duplicates. @@ -135,39 +86,50 @@ find_nondominated_3d_impl_sorted(const double ** restrict rows, size_t size, } early_end: - free(tnodes); + treap2d_free(tree); return new_size; } - /** - Returns NULL if point is dominated by a different point in the tree. + Return true if point is dominated by a different point in the tree. */ -static avl_node_t * -dominated_by_tree_3d(const avl_tree_t * restrict tree, - const double * restrict point) +static inline bool +dominated_by_tree_3d(Treap2D *tree, const double *point) { - avl_node_t * nodeaux; - int res = avl_search_closest(tree, point, &nodeaux); - if (res >= 0) { // nodeaux goes before point - const double * restrict prev = nodeaux->item; - DEBUG2(printf_point("res > 0: prev = [ ", prev, 3, " ], ")); - DEBUG2(printf_point("point = [ ", point, 3, " ]\n")); - assert(prev[0] != INFINITY); - assert(prev[0] <= point[0]); - assert(prev[1] <= point[1]); - return prev[2] <= point[2] ? NULL : nodeaux->next; - } else if (nodeaux->prev) { // nodeaux goes after point, so move to the next one. - const double * restrict prev = nodeaux->prev->item; - DEBUG2(printf_point("res <= 0: prev = [ ", prev, 3, " ], ")); - DEBUG2(printf_point("point = [ ", point, 3, " ]\n")); - assert(prev[0] != INFINITY); - assert(prev[0] <= point[0]); - assert(prev[1] <= point[1]); - return prev[2] <= point[2] ? NULL : nodeaux; + DEBUG1(treap2d_validate_tree(tree)); + point++; // Ignore dim [0] + // Find the rightmost point with prev->key <= point[0]. + TreapNode *prev_node = treap2d_find_le(tree, point[0]); + if (prev_node != NULL && treap2d_node_get_y(prev_node) <= point[1]) { + assert(treap2d_node_get_x(prev_node) <= point[0]); + return true; // p is dominated } - assert(nodeaux != NULL); - return nodeaux; // point is not dominated + return false; // p is not dominated +} + +// FIXME: How to merge this function and the one above? +static inline bool +insert_if_not_dominated_by_tree_3d(Treap2D *tree, const double *point) +{ + DEBUG1(treap2d_validate_tree(tree)); + point++; // Ignore dim [0] + // Find the rightmost point with prev->key <= point[0]. + TreapNode **prev_link = treap2d_find_le_link(tree, point[0]); + if (prev_link != NULL) { + const double prev[2] = { treap2d_node_get_x(*prev_link), treap2d_node_get_y(*prev_link) }; + if (prev[1] <= point[1]) { + assert(prev[0] <= point[0]); + return false; // p is dominated + } else if (prev[0] == point[0]) { + // p dominates prev. Remove the already found node. + treap2d_erase_at(prev_link); + DEBUG1(treap2d_validate_tree(tree)); + } + } + // FIXME: This will call split_lt, but we just called find_le above, so can + // we avoid one of them? + treap2d_insert(tree, point[0], point[1]); + return true; // p was inserted. } /** @@ -192,17 +154,10 @@ kung_merge_dim3(const double ** restrict r, size_t r_size, DEBUG1(for (size_t j = 0; j < k; j++) assert(s[j][0] < r0)); DEBUG1(for (size_t j = k; j < s_size; j++) assert(r0 <= s[j][0])); - avl_tree_t tree; - avl_init_tree(&tree, qsort_cmp_pdouble_asc_y_asc_z); - // FIXME: Use a workspace to allocate this once and re-alloc only if a larger number is needed. - avl_node_t * tnodes = malloc((r_size + 1) * sizeof(*tnodes)); - avl_node_t * node = tnodes; - node->item = r[0]; - avl_insert_top(&tree, node); - - const double sentinel[] = { INFINITY, INFINITY, -INFINITY}; - (++node)->item = sentinel; - avl_insert_after(&tree, node - 1, node); + Treap2D * tree = treap2d_new(r_size); + assert(tree != NULL); + treap2d_init_with_single_node(tree, r[0][1], r[0][2]); + DEBUG2(printf_point("insert in tree: r=[ ", r[0], 3, " ]\n")); size_t i = 1, new_size = s_size; do { @@ -212,39 +167,21 @@ kung_merge_dim3(const double ** restrict r, size_t r_size, while (i < r_size) { // Add to the tree all points in R that could dominate v. const double * restrict u = r[i]; - if (u[0] > v[0]) { + if (u[0] > v[0]) break; - } - avl_node_t * nodeaux = dominated_by_tree_3d(&tree, u); - if (nodeaux != NULL) { // u is NOT dominated by a point in the tree. - const double * restrict point = nodeaux->item; - assert(u[1] <= point[1]); - // Delete everything in the tree that is dominated by u. - while (u[2] <= point[2]) { - // printf("delete point: "); print_point(point); printf("\n"); - assert(u[1] <= point[1]); - nodeaux = nodeaux->next; - point = nodeaux->item; - /* FIXME: A possible speed up is to delete without - rebalancing the tree because avl_insert_before() will - rebalance, but we need to know which is the highest node - that needs rebalancing. */ - avl_unlink_node(&tree, nodeaux->prev); - } - // printf("insert before point: "); print_point(point); printf("\n"); - (++node)->item = u; - avl_insert_before(&tree, nodeaux, node); + if (insert_if_not_dominated_by_tree_3d(tree, u)) { + // u is NOT dominated by a point in the tree. + DEBUG2(printf_point("!dominated_by_tree_3d: u=[ ", u, 3, " ]\n")); } i++; } - if (dominated_by_tree_3d(&tree, v) == NULL) { + if (dominated_by_tree_3d(tree, v)) { DEBUG2(printf_point("dominated_by_tree_3d: v=[ ", v, 3, " ]\n")); s[k] = NULL; // dominated new_size--; } k++; } while (k < s_size); - - free(tnodes); + treap2d_free(tree); return new_size; } diff --git a/c/pareto.c b/c/pareto.c index 42e7dbf67..b1e1a27f6 100644 --- a/c/pareto.c +++ b/c/pareto.c @@ -4,18 +4,7 @@ #include "io.h" #endif -typedef const double avl_item_t; -typedef struct avl_node_t { - struct avl_node_t *next; - struct avl_node_t *prev; - struct avl_node_t *parent; - struct avl_node_t *left; - struct avl_node_t *right; - avl_item_t *item; - unsigned char depth; -} avl_node_t; - -#include "avl_tiny.h" +#include "treap2d_fixed_size.h" /** Nondominated sorting in 3D in O(k * n log n), where k is the number of fronts. @@ -33,21 +22,14 @@ pareto_rank_3d(int * restrict rank, const double * restrict points, size_t size) const bool keep_weakly = true; const double ** p = generate_row_pointers_asc_rev_3d(points, size); - avl_tree_t tree; - avl_init_tree(&tree, qsort_cmp_pdouble_asc_x_nonzero); - avl_node_t * tnodes = malloc((size+1) * sizeof(*tnodes)); - const double sentinel[] = { INFINITY, -INFINITY }; - tnodes->item = sentinel; + Treap2D * tree = treap2d_new(size); + assert(tree != NULL); int front = 0; while (true) { ASSUME(size >= 2); const double * restrict pk = p[0]; - avl_node_t * node = tnodes + 1; - node->item = pk; - avl_insert_top(&tree, node); - // Insert sentinel. - avl_insert_after(&tree, node, tnodes); + treap2d_init_with_single_node(tree, pk[0], pk[1]); // In this context, size means "no dominated solution found". size_t n_nondom = size, j = 1; @@ -56,41 +38,10 @@ pareto_rank_3d(int * restrict rank, const double * restrict points, size_t size) const double * restrict pj = p[j]; bool dominated; if (pk[0] > pj[0] || pk[1] > pj[1]) { - avl_node_t * nodeaux; - int res = avl_search_closest(&tree, pj, &nodeaux); - assert(res != 0); - if (res > 0) { // nodeaux goes before pj - const double * restrict prev = nodeaux->item; - assert(prev[0] != sentinel[0]); - assert(prev[0] <= pj[0]); - dominated = prev[1] <= pj[1]; - 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]; - } else { - dominated = false; - } - - if (!dominated) { // pj is NOT dominated by a point in the tree. - const double * restrict point = nodeaux->item; - assert(pj[0] <= point[0]); - // Delete everything in the tree that is dominated by pj. - while (pj[1] <= point[1]) { - assert(pj[0] <= point[0]); - nodeaux = nodeaux->next; - point = nodeaux->item; - /* FIXME: A possible speed up is to delete without - rebalancing the tree because avl_insert_before() - will rebalance, but we need to know which is the - highest node that needs rebalancing. */ - avl_unlink_node(&tree, nodeaux->prev); - } - (++node)->item = pj; - avl_insert_before(&tree, nodeaux, node); - } + TreapNode *pred = treap2d_find_le(tree, pj[0]); + dominated = (pred != NULL && treap2d_node_get_y(pred) <= pj[1]); + if (!dominated) + treap2d_insert(tree, pj[0], pj[1]); } else { // Handle duplicates and points that are dominated by the immediate // previous one. @@ -126,7 +77,7 @@ pareto_rank_3d(int * restrict rank, const double * restrict points, size_t size) // we can stop. size -= n_nondom; if (size <= 1) { - free(tnodes); + treap2d_free(tree); free(p); return; } @@ -147,7 +98,6 @@ pareto_rank_3d(int * restrict rank, const double * restrict points, size_t size) k, p[k][0], p[k][1], p[k][2]); front++; - avl_clear_tree(&tree); } } diff --git a/c/solutions_list.h b/c/solutions_list.h new file mode 100644 index 000000000..1c819d9ab --- /dev/null +++ b/c/solutions_list.h @@ -0,0 +1,149 @@ +#ifndef SOLUTIONS_LIST_H_ +#define SOLUTIONS_LIST_H_ + +#include +#include // memcpy +#include + +typedef struct VoidList { + uint32_t n; + uint32_t cap; + const void ** list; +} VoidList; + +/** + A duplicate-aware leaf bucket stores one unique objective vector plus all + payloads that share it, so duplicate-heavy leaves are checked once per + distinct objective vector rather than once per stored solution. +*/ +typedef struct SolutionsList { + const double * z; // point coordinates (objective vector). + // FIXME: Compile with -fms-extensions to avoid duplication and use: + // struct VoidList; + VoidList x; +} SolutionsList; + +static inline bool +void_list_init(VoidList * restrict item, const void * restrict x) +{ + item->cap = 1; + item->n = 1; + if (x == NULL) { + item->list = NULL; + } else { + item->list = (const void **) malloc(item->cap * sizeof(*item->list)); + if (item->list == NULL) + return false; + item->list[0] = x; + } + return true; +} + +static inline bool +solutions_list_init(SolutionsList * restrict item, const double * restrict z, const void * restrict x, dimension_t dim) +{ + assert(dim >= 2); + // FIXME: We should have a pool of z vectors or make z a flexible array. + double * new_z = malloc(dim * sizeof(*new_z)); + if (new_z == NULL) + return false; + item->z = memcpy(new_z, z, dim * sizeof(*z)); + return void_list_init(&item->x, x); +} + +static inline void +void_list_free(VoidList * item) +{ + item->n = 0; + // item->n may be > 0 if we stored repeated z values but no x values. + if (!item->list) + return; + free(item->list); + item->list = NULL; +} + +static inline void +solutions_list_free(SolutionsList * item) +{ + free((void *)item->z); + item->z = NULL; + void_list_free(&item->x); +} + +/** + Add a duplicated point. +*/ +static inline bool +void_list_append(VoidList * restrict item, const void * restrict x) +{ + if (x == NULL) { + item->n++; + return true; + } + + if (item->n == item->cap) { + assert(item->cap < UINT32_MAX / 2); + uint32_t new_cap = item->cap * 2; + const void ** new_list = (const void **) realloc(item->list, new_cap * sizeof(*new_list)); + if (new_list == NULL) + return false; + + item->list = new_list; + item->cap = new_cap; + } + item->list[item->n++] = x; + return true; +} + +/** + Add a duplicated point. +*/ +static inline bool +solutions_list_append(SolutionsList * restrict item, const void * restrict x) +{ + assert(item->z != NULL); + return void_list_append(&item->x, x); +} + + +static inline void +void_list_get_x_values(const VoidList * item, const void *** x_list) +{ + for (unsigned i = 0; i < item->n; i++, (*x_list)++) { + assert(item->list[i] != NULL); // Storing NULL values is a bug. + **x_list = item->list[i]; + } +} + +static inline void +solutions_list_get_x_values(const SolutionsList * item, const void *** x_list) +{ + void_list_get_x_values(&item->x, x_list); +} + +static inline void +solutions_list_get_contents(const SolutionsList * item, + double ** z_list, const void *** x_list, + dimension_t dim) +{ + for (unsigned i = 0; i < item->x.n; i++, *z_list += dim) { + // Duplicate z for each x + memcpy(*z_list, item->z, dim * sizeof **z_list); + } + if (x_list != NULL) + void_list_get_x_values(&item->x, x_list); +} + + +/** + Shallow copy. +*/ +static inline void +solutions_list_copy(SolutionsList * restrict dst, const SolutionsList * restrict src) +{ + memcpy(dst, src, sizeof *src); +} + + + +#endif // SOLUTIONS_LIST_H_ diff --git a/c/sort.h b/c/sort.h index c819c57a9..a16d297fe 100644 --- a/c/sort.h +++ b/c/sort.h @@ -91,6 +91,21 @@ all_equal_double(const double * restrict a, const double * restrict b, dimension return (bool) a_eq_b; } +/** + Return true if a <= b but a != b, false otherwise. +*/ +static inline bool +dominates_2d(const double * restrict a, const double * restrict b) +{ + return a[0] <= b[0] && a[1] <= b[1] && (a[0] != b[0] || a[1] != b[1]); +} + +static inline bool +dominates(const double * restrict a, const double * restrict b, dimension_t dim) +{ + return weakly_dominates(a, b, dim) && !all_equal_double(a, b, dim); +} + // ---------- Comparison functions (e.g, qsort). Return 'int' ---------------- // General type for comparison functions used in qsort() and qsort_r(). diff --git a/c/treap2d.h b/c/treap2d.h new file mode 100644 index 000000000..4fe65c8e2 --- /dev/null +++ b/c/treap2d.h @@ -0,0 +1,886 @@ +/***************************************************************************** + + treap.h: Randomised binary search tree (Treap) for a 2D mutually nondominated + frontier. + + --------------------------------------------------------------------- + Copyright (C) 2026 + Manuel Lopez-Ibanez + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + --------------------------------------------------------------------- + + [1] R. Seidel and C. R. Aragon. Randomized search trees. Algorithmica, + 16:464-497, 1996 + +******************************************************************************/ +#ifndef TREAP_H_ +#define TREAP_H_ + +#include "common.h" +#include "sort.h" + +/** + + Randomized binary search tree for a 2-D mutually non-dominated + minimization frontier. + + Frontier invariant: + + x increases strictly from left to right + y decreases strictly from left to right + + Consequently, for a new point (x,y), the old frontier points + dominated by the new point form a contiguous range in key order. + This allows the whole dominated range to be detached with + split/merge operations instead of deleting its nodes individually. + + Nodes are externally owned. The treap does not allocate or free them. + + Expected complexity: + + search O(log n) + split O(log n) + merge O(log n) + insert O(log n) + insert+displace O(log n) + + Worst-case complexity of a treap operation is O(n), as for any + randomized BST. +*/ +#include +#include +#include +#include +#include + +struct TreapNode { + double x, y; + uint64_t priority; // Parents should have a higher priority than children. + struct TreapNode *left; + struct TreapNode *right; +#ifndef TREAP_ITEM_IS_EMPTY + TreapItem item; +#endif +}; +typedef struct TreapNode TreapNode; + +struct Treap { + TreapNode *root; + uint64_t rng_state; +}; +typedef struct Treap Treap; + + +/** + PRNG is xorshift64*, which is not cryptographically safe but it is fast. +*/ +static inline uint64_t +treap_random(Treap *tree) +{ + uint64_t x = tree->rng_state; + assert(x != 0); // State must never be zero. + + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + + tree->rng_state = x; + return x * UINT64_C(0x2545F4914F6CDD1D); +} + + +/** + Treap Seed generation + + Use the address of the Treap object as process-dependent entropy, then + thoroughly mix it with SplitMix64's mixing function. + + This is not intended as cryptographic randomness. Its purpose is simply + to give separately allocated treaps normally different PRNG streams, + without requiring global mutable state or an OS RNG call. +*/ +static inline uint64_t +treap_seed(const Treap *tree) +{ + // Cast to uintptr_t to use the pointer value as an integer. + uint64_t x = (uint64_t)(uintptr_t)tree; + /* Mix the address. The addition is useful even though a fixed constant + would ultimately be mixed anyway. */ + x += UINT64_C(0x9e3779b97f4a7c15); + + x = (x ^ (x >> 30)) * UINT64_C(0xbf58476d1ce4e5b9); + x = (x ^ (x >> 27)) * UINT64_C(0x94d049bb133111eb); + x ^= x >> 31; + // xorshift64* cannot leave the zero state. + return x ? x : UINT64_C(0x9e3779b97f4a7c15); +} + + +/** + (Re)-Initializes the tree structure. Nothing is free()d. + + Complexity: O(1) +*/ +static inline void +treap_init(Treap *tree) +{ + tree->root = NULL; + tree->rng_state = treap_seed(tree); +} + +static inline void +treap_init_with_single_node(Treap *tree, TreapNode * node) +{ + assert(node->left == NULL); + assert(node->right == NULL); + treap_init(tree); + tree->root = node; + node->priority = treap_random(tree); +} + +static inline void +treap_node_init_empty(TreapNode *node, double x, double y) +{ + assert(node != NULL); + node->x = x; + node->y = y; + node->left = NULL; + node->right = NULL; +} + +#ifndef TREAP_ITEM_IS_EMPTY +/** + Initialize a node before inserting it. + + The node does not receive a priority until insertion because the + priority should come from the destination treap's RNG. +*/ +static inline void +treap_node_init(TreapNode *node, double x, double y, TreapItem item) +{ + treap_node_init_empty(node, x, y); + node->item = item; +} + +static inline void +treap_init_with_singleton(Treap *tree, TreapNode * node, double x, double y, TreapItem item) +{ + treap_node_init(node, x, y, item); + treap_init_with_single_node(tree, node); +} +#else +static inline void +treap_node_init(TreapNode *node, double x, double y) +{ + treap_node_init_empty(node, x, y); +} + +static inline void +treap_init_with_singleton(Treap *tree, TreapNode * node, double x, double y) +{ + treap_node_init(node, x, y); + treap_init_with_single_node(tree, node); +} +#endif + + + + +static inline void +treap_split_lt_with_max(TreapNode *root, double key, + TreapNode **left, TreapNode **right, + TreapNode **left_max) +{ + if (root == NULL) { + *left = NULL; + *right = NULL; + if (left_max) + *left_max = NULL; + return; + } + + if (root->x < key) { + *left = root; + if (left_max) { + TreapNode *sub_left_max; + treap_split_lt_with_max(root->right, key, &root->right, right, + &sub_left_max); + /* If anything remained in root->right, its maximum is the maximum of + the whole left result. Otherwise root is the max. */ + *left_max = sub_left_max != NULL ? sub_left_max : root; + } else { + treap_split_lt_with_max(root->right, key, &root->right, right, NULL); + } + } else { + *right = root; + treap_split_lt_with_max(root->left, key, left, &root->left, left_max); + } +} + +/** + Split by x. + + On return: + + *left contains every node with x < key + *right contains every node with x >= key + + The frontier and treap invariants are preserved. + + Expected complexity: O(log n). +*/ +static inline void +treap_split_lt(TreapNode *root, double key, + TreapNode **left, TreapNode **right) +{ + treap_split_lt_with_max(root, key, left, right, NULL); +} + + +/** + Split by key. + + On return: + + *left contains every node with x <= key + *right contains every node with x > key + + The frontier and treap invariants are preserved. + + Expected complexity: O(log n). +*/ +static inline void +treap_split_le(TreapNode *root, double key, + TreapNode **left, TreapNode **right) +{ + if (root == NULL) { + *left = NULL; + *right = NULL; + return; + } + + if (root->x <= key) { + // root belongs to the left result. + *left = root; + treap_split_le(root->right, key, &root->right, right); + } else { + // root belongs to the right result. + *right = root; + treap_split_le(root->left, key, left, &root->left); + } +} + +/** + Split a frontier treap by z[1]. + + Since z[1] strictly decreases as z[0] increases: + + a: z[1] >= value + b: z[1] < value + + The nodes in `a` therefore form a prefix in key order. + + Complexity: O(log n) assuming a balanced tree. +*/ +static inline void +treap_split_z1_ge_with_min(TreapNode *root, double y, + TreapNode **high_y, TreapNode **low_y, + TreapNode **low_y_min) +{ + if (root == NULL) { + *high_y = NULL; + *low_y = NULL; + if (low_y_min) + *low_y_min = NULL; + return; + } + + if (root->y >= y) { + *high_y = root; + treap_split_z1_ge_with_min(root->right, y, &root->right, low_y, low_y_min); + } else { + *low_y = root; + if (low_y_min) { + TreapNode *sub_low_min; + treap_split_z1_ge_with_min(root->left, y, high_y, &root->left, &sub_low_min); + /* If root has a left part remaining, its minimum is the minimum of + that part. Otherwise root is the minimum. */ + if (low_y_min) + *low_y_min = sub_low_min != NULL ? sub_low_min : root; + } else { + treap_split_z1_ge_with_min(root->left, y, high_y, &root->left, NULL); + } + } +} + +/** + Split a frontier by y. + + PRECONDITION: + + root is a valid mutually non-dominated frontier. + + Since y strictly decreases as x increases: + + *high_y contains every node with node->y >= y + *low_y contains every node with node->y < y + + Therefore *high_y is a prefix in x order. + + This is the operation used to detach all nodes dominated by a new + frontier point. + + Expected complexity: O(log n). +*/ +static inline void +treap_split_z1_ge(TreapNode *root, double y, + TreapNode **high_y, TreapNode **low_y) +{ + treap_split_z1_ge_with_min(root, y, high_y, low_y, NULL); +} + +static TreapNode * +treap_rotate_left(TreapNode *node) +{ + TreapNode *right = node->right; + node->right = right->left; + right->left = node; + return right; +} + +static TreapNode * +treap_rotate_right(TreapNode *node) +{ + TreapNode *left = node->left; + node->left = left->right; + left->right = node; + return left; +} + +static TreapNode * +treap_insert_node_and_rotate(TreapNode * root, TreapNode *node) +{ + if (!root) + return node; + + if (node->x < root->x) { + root->left = treap_insert_node_and_rotate(root->left, node); + assert(root->left != NULL); + if (root->left->priority > root->priority) + root = treap_rotate_right(root); + } else { + root->right = treap_insert_node_and_rotate(root->right, node); + assert(root->right != NULL); + if (root->right->priority > root->priority) + root = treap_rotate_left(root); + } + return root; +} + +/** + Merge two treaps. + + PRECONDITION: + + Every x in left is less than every x in right. + + Returns the merged treap. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_merge(TreapNode *left, TreapNode *right) +{ + if (!left) + return right; + + if (!right) + return left; + + if (left->priority > right->priority) { + left->right = treap_merge(left->right, right); + return left; + } + + right->left = treap_merge(left, right->left); + return right; +} + +static inline void +treap_validate_tree(TreapNode * node) +{ + if (!node) + return; + + if (node->left) { + assert(node->left->x < node->x); + // Ideally, there should be no equal priorities, but it does happen. + assert(node->priority >= node->left->priority); + treap_validate_tree(node->left); + } + if (node->right) { + assert(node->right->x > node->x); + // Ideally, there should be no equal priorities, but it does happen. + assert(node->priority >= node->right->priority); + treap_validate_tree(node->right); + } +} + +// O(log n) +static inline const TreapNode * +treap_minimum(const TreapNode * root) +{ + assert(root != NULL); + while (root->left != NULL) + root = root->left; + + return root; +} + +// O(log n) +static inline const TreapNode * +treap_maximum(const TreapNode *root) +{ + assert(root != NULL); + + while (root->right != NULL) + root = root->right; + + return root; +} + +/** + Find the node with the largest x such that x <= key. + + Returns NULL if no such node exists. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_find_le(const Treap *tree, double key) +{ + TreapNode *node = tree->root, *best = NULL; + + while (node != NULL) { + if (node->x <= key) { + best = node; + node = node->right; + } else { + node = node->left; + } + } + return best; +} + +static inline TreapNode ** +treap_find_le_link(Treap *tree, double key) +{ + TreapNode **link = &tree->root, **best = NULL; + + while (*link != NULL) { + TreapNode *node = *link; + + if (node->x <= key) { + best = link; + link = &node->right; + } else { + link = &node->left; + } + } + + return best; +} + +/** + Find the node with the largest x such that x < key. + + Returns NULL if no such node exists. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_find_lt(const Treap *tree, double key) +{ + TreapNode *node = tree->root, *best = NULL; + + while (node != NULL) { + if (node->x < key) { + best = node; + node = node->right; + } else { + node = node->left; + } + } + return best; +} +/** + Find the node with the smallest x such that x >= key. + + Returns NULL if no such node exists. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_find_ge(const Treap *tree, double key) +{ + TreapNode *node = tree->root, *best = NULL; + + while (node != NULL) { + if (node->x >= key) { + best = node; + node = node->left; + } else { + node = node->right; + } + } + return best; +} + +/** + Find the node with the smallest x such that x > key. + + Returns NULL if no such node exists. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_find_gt(const Treap *tree, double key) +{ + TreapNode *node = tree->root, *best = NULL; + + while (node != NULL) { + if (node->x > key) { + best = node; + node = node->left; + } else { + node = node->right; + } + } + return best; +} + + +/** + Find a node with exactly x == key. + + Returns NULL if absent. + + Expected complexity: O(log n). + */ +static inline TreapNode * +treap_find(const Treap *tree, double key) +{ + TreapNode *node = tree->root; + + while (node) { + if (key < node->x) + node = node->left; + else if (key > node->x) + node = node->right; + else + return node; + } + return NULL; +} + +static inline bool +treap_dominates(const Treap * tree, const double * z) +{ + TreapNode * node = treap_find_le(tree, z[0]); + return node && dominates_2d((double[2]){node->x, node->y}, z); +} + +/** + Insert a single node. + + PRECONDITIONS: + + - node is initialized with treap_node_init(). + - node->x does not already exist in the treap. + - node->left == NULL. + - node->right == NULL. + + This function does not check whether the resulting set is a valid + Pareto frontier. Use treap_insert_and_displace() when inserting + a new non-dominated point into an existing frontier. +*/ +static inline void +treap_insert_node(Treap *tree, TreapNode *node) +{ + assert(tree != NULL); + assert(node != NULL); + assert(node->left == NULL); + assert(node->right == NULL); + // node->x is not already in the tree. + assert(treap_find(tree, node->x) == NULL); + DEBUG1(treap_validate_tree(tree->root)); + + node->priority = treap_random(tree); + +#if 0 + TreapNode *left; + TreapNode *right; + /* Since keys are guaranteed unique by the caller, there is no need to + search for an existing x before doing the split. */ + treap_split_lt(tree->root, node->x, &left, &right); + tree->root = treap_merge(treap_merge(left, node), right); +#else + tree->root = treap_insert_node_and_rotate(tree->root, node); +#endif + DEBUG1(treap_validate_tree(tree->root)); +} + +static inline void +treap_node_free(TreapNode * node) +{ + if (node == NULL) + return; + + treap_node_free(node->left); + node->left = NULL; + treap_node_free(node->right); + node->right = NULL; +#ifndef TREAP_ITEM_IS_EMPTY + treap_item_free(node->item); +#endif + free(node); +} + +static inline void +treap_free_nodes(Treap * t) +{ + treap_node_free(t->root); + t->root = NULL; +} + +static inline TreapNode * +treap_erase_at(TreapNode **link) +{ + assert(link != NULL); + assert(*link != NULL); + + TreapNode *node = *link; + *link = treap_merge(node->left, node->right); + node->left = NULL; + node->right = NULL; + return node; +} + +/** + Inserts node into tree and detaches all nodes dominated by node. + + On return: + + *prev = node's predecessor in the resulting tree + *next = node's successor in the resulting tree + + Returns the root of a separate treap containing the displaced nodes. + + Expected complexity: O(log n). + + FIXME: This function performs 4 tree traversals. It may be possible to do + it faster while rotating nodes as needed. +*/ +static inline TreapNode * +treap_insert_and_displace_get_bounds(Treap *tree, TreapNode *node, + TreapNode **prev, TreapNode **next) +{ + assert(tree != NULL); + assert(node != NULL); + assert(node->left == NULL); + assert(node->right == NULL); + + const double x = node->x; + const double y = node->y; + TreapNode *left, *right; + if (prev) { + TreapNode *left_max; + treap_split_lt_with_max(tree->root, x, &left, &right, &left_max); + *prev = left_max; + } else { + treap_split_lt(tree->root, x, &left, &right); + } + + /* On the frontier, z[1] decreases with increasing x. Therefore the + nodes in right with z[1] >= y are exactly the nodes dominated by + node, and they form a contiguous prefix. */ + TreapNode *displaced, *keep; + if (next) { + TreapNode *keep_min; + treap_split_z1_ge_with_min(right, y, &displaced, &keep, &keep_min); + *next = keep_min; + } else { + treap_split_z1_ge(right, y, &displaced, &keep); + } + + // Insert the new node between left and keep. + node->priority = treap_random(tree); + node->left = NULL; + node->right = NULL; + // key ordering is: left < node < keep + tree->root = treap_merge(treap_merge(left, node), keep); + DEBUG1(treap_validate_tree(tree->root)); + return displaced; +} + +/** + Insert a new point into a mutually non-dominated frontier and detach all + existing nodes dominated by it. + + PRECONDITIONS: + + - tree is a valid mutually non-dominated frontier. + - node is not already in tree. + - node->x is unique in tree. + - no node in tree weakly dominates node. + - node->left == NULL. + - node->right == NULL. + + POSTCONDITIONS: + + - node belongs to tree. + - the returned root is a separate treap containing exactly the nodes + removed from tree. + - every returned node is dominated by node. + - no node other than those dominated by node is removed. + + Expected complexity: O(log n). +*/ +static inline TreapNode * +treap_insert_and_displace(Treap *tree, TreapNode *node) +{ + return treap_insert_and_displace_get_bounds(tree, node, NULL, NULL); +} + +/** + Insert the mutually-non-dominated frontier `incoming` into `existing`. + + The incoming nodes are allowed to have priorities unrelated to + `existing`; every incoming node receives a fresh priority. + + Returns the resulting frontier and places all displaced existing + nodes in *displaced. + + Preconditions: + + - incoming is a valid frontier: + z[0] strictly increasing + z[1] strictly decreasing + + - existing is a valid frontier with the same property + + - keys are unique within each tree + + - no node in existing weakly dominates any node in incoming + + Postconditions: + + - result is a valid frontier + - *displaced is a valid frontier + - every displaced node is from existing + - every existing node dominated by an incoming node is displaced + - no incoming node is displaced +*/ +static inline TreapNode * +treap_frontier_union(Treap *t, TreapNode *incoming, TreapNode *existing, + TreapNode **displaced, bool *all_displaced) +{ + if (!incoming) { + *displaced = NULL; + *all_displaced = (existing == NULL); + return existing; + } + + if (!existing) { + *displaced = NULL; + *all_displaced = true; + return incoming; + } + + /* Take the incoming root as the pivot. It is guaranteed not to be + dominated by any node in existing. */ + const double x = incoming->x; + const double y = incoming->y; + + TreapNode *existing_left; // existing->x < x + TreapNode *existing_ge; // existing->x >= x + treap_split_lt(existing, x, &existing_left, &existing_ge); + + TreapNode *equal; // existing_ge->x == x + TreapNode *existing_right; // existing_ge->x > x + treap_split_le(existing_ge, x, &equal, &existing_right); + + /* Since y decreases as x increases, the nodes in existing_right + dominated by (x,y) form its prefix: existing_right->y >= y */ + TreapNode *dominated; + TreapNode *existing_keep; + treap_split_z1_ge(existing_right, y, &dominated, &existing_keep); + + // Recursively process the two independent key ranges. + TreapNode *displaced_left, *displaced_right; + bool all_displaced_left, all_displaced_right; + TreapNode *left = treap_frontier_union(t, incoming->left, existing_left, + &displaced_left, &all_displaced_left); + TreapNode *right = treap_frontier_union(t, incoming->right, existing_keep, + &displaced_right, &all_displaced_right); + + *all_displaced = all_displaced_left && all_displaced_right; + + /* The displaced nodes are also already separated by key: + + displaced_left < equal < dominated < displaced_right + + `equal` is either NULL or a singleton. + */ + // FIXME: Can we move this up? + *displaced = treap_merge(treap_merge(displaced_left, equal), + treap_merge(dominated, displaced_right)); + + // The incoming node is now a singleton. + incoming->priority = treap_random(t); + incoming->left = NULL; + incoming->right = NULL; + + /* Reconstruct the incoming result. + + The two recursive results are ordered: + + left < incoming < right + + Ordinary treap merges restore the heap property regardless of + the incoming node's new priority. + */ + return treap_merge(treap_merge(left, incoming), right); +} + + +static inline void +treap_node_get_unique_vectors(const TreapNode *node, const double *** z_list) +{ + if (!node) + return; + + treap_node_get_unique_vectors(node->left, z_list); + **z_list = &node->x; + (*z_list)++; + treap_node_get_unique_vectors(node->right, z_list); +} + +/** + z_list must be an array of pointers to 'const double' of length treap_unique_size(t). + + That is, z_list does not contain the actual data but just pointers into the + Treap. The data is owned by the Treap, so it should not be modified using + these pointers. +*/ +static inline void +treap_get_unique_vectors(const Treap *tree, const double ** z_list) +{ + treap_node_get_unique_vectors(tree->root, &z_list); +} + +#endif /* TREAP_H_ */ diff --git a/c/treap2d_dlnode.h b/c/treap2d_dlnode.h new file mode 100644 index 000000000..b596a0b3e --- /dev/null +++ b/c/treap2d_dlnode.h @@ -0,0 +1,74 @@ +#ifndef TREAP_2D_DLNODE_H +#define TREAP_2D_DLNODE_H + +#include +#include "gcc_attribs.h" + +typedef dlnode_t * TreapItem; +static inline void treap_item_free(TreapItem item _attr_maybe_unused) { return; } + +#include "treap2d.h" + +typedef struct Treap2D { + Treap treap; + TreapNode * last_node; + TreapNode node_pool[]; // Flexible array +} Treap2D; + +static Treap2D * +treap2d_new(size_t n) +{ + Treap2D * tree = malloc(sizeof(Treap2D) + n * sizeof(TreapNode)); + tree->last_node = tree->node_pool; + return tree; +} + +static void +treap2d_free(Treap2D * tree) +{ + free(tree); +} + +static inline void +treap2d_init_node(TreapNode * restrict node, dlnode_t * restrict p) +{ + // The treap is sorted in ascending p->x[1] and descending p->x[0]. + treap_node_init(node, p->x[1], p->x[0], p); +} + +static inline dlnode_t * +treap2d_get_dlnode(TreapNode * node) +{ + return node->item; +} + +static void +treap2d_init_with_sentinels(Treap2D *tree, dlnode_t * restrict p, dlnode_t * restrict s2, dlnode_t * restrict s1) +{ + TreapNode * node = tree->last_node; + treap2d_init_node(node, p); + treap2d_init_node(node + 1, s2); + (node+1)->priority = 0; // Push to the bottom. + treap2d_init_node(node + 2, s1); + (node+2)->priority = 0; // Push to the bottom. + treap_init_with_single_node(&tree->treap, node); + node->left = node + 1; + node->right = node + 2; + tree->last_node += 3; +} + +static TreapNode * +treap2d_find_le(const Treap2D *tree, double value) +{ + return treap_find_le(&tree->treap, value); +} + +static inline void +treap2d_insert_and_displace_get_bounds(Treap2D *tree, dlnode_t * restrict p, TreapNode **prev, TreapNode **next) +{ + treap2d_init_node(tree->last_node, p); + (void)treap_insert_and_displace_get_bounds(&tree->treap, tree->last_node, prev, next); + tree->last_node++; +} + +#endif // TREAP_2D_DLNODE_H diff --git a/c/treap2d_fixed_size.c b/c/treap2d_fixed_size.c new file mode 100644 index 000000000..54080dec6 --- /dev/null +++ b/c/treap2d_fixed_size.c @@ -0,0 +1,77 @@ +#include "treap2d_fixed_size.h" + +#define TREAP_ITEM_IS_EMPTY +#include "treap2d.h" + +typedef struct Treap2D { + Treap treap; + TreapNode * last_node; + TreapNode node_pool[]; // Flexible array +} Treap2D; + + +Treap2D * +treap2d_new(size_t n) +{ + Treap2D * tree = malloc(sizeof(Treap2D) + n * sizeof(TreapNode)); + tree->last_node = tree->node_pool; + return tree; +} + +void +treap2d_free(Treap2D * tree) +{ + free(tree); +} + +void +treap2d_init_with_single_node(Treap2D *tree, double x, double y) +{ + treap_node_init(tree->last_node, x, y); + treap_init_with_single_node(&tree->treap, tree->last_node); + tree->last_node++; +} + +TreapNode * +treap2d_find_le(const Treap2D *tree, double value) +{ + return treap_find_le(&tree->treap, value); +} + +TreapNode ** +treap2d_find_le_link(Treap2D *tree, double value) +{ + return treap_find_le_link(&tree->treap, value); +} + +double +treap2d_node_get_x(const TreapNode *node) +{ + return node->x; +} + +double +treap2d_node_get_y(const TreapNode *node) +{ + return node->y; +} + +void +treap2d_validate_tree(Treap2D * tree) +{ + treap_validate_tree(tree->treap.root); +} + +TreapNode * +treap2d_erase_at(TreapNode **link) +{ + return treap_erase_at(link); +} + +void +treap2d_insert(Treap2D *tree, double x, double y) +{ + treap_node_init(tree->last_node, x, y); + treap_insert_and_displace(&tree->treap, tree->last_node); + tree->last_node++; +} diff --git a/c/treap2d_fixed_size.h b/c/treap2d_fixed_size.h new file mode 100644 index 000000000..2f9b66877 --- /dev/null +++ b/c/treap2d_fixed_size.h @@ -0,0 +1,20 @@ +#ifndef TREAP_2D_FIXED_SIZE_H_ +#define TREAP_2D_FIXED_SIZE_H_ + +#include + +typedef struct Treap2D Treap2D; +typedef struct TreapNode TreapNode; + +Treap2D * treap2d_new(size_t n); +void treap2d_free(Treap2D * tree); +void treap2d_init_with_single_node(Treap2D *tree, double x, double y); +TreapNode * treap2d_find_le(const Treap2D *tree, double value); +TreapNode ** treap2d_find_le_link(Treap2D *tree, double value); +double treap2d_node_get_x(const TreapNode *node); +double treap2d_node_get_y(const TreapNode *node); +void treap2d_validate_tree(Treap2D * tree); +TreapNode * treap2d_erase_at(TreapNode **link); +void treap2d_insert(Treap2D *tree, double x, double y); + +#endif // TREAP_2D_FIXED_SIZE_H_ diff --git a/c/treap_archive.c b/c/treap_archive.c new file mode 100644 index 000000000..25a3ac141 --- /dev/null +++ b/c/treap_archive.c @@ -0,0 +1,390 @@ +/************************************************************************* + + treap_archive.c: 2D Archive based on Randomised Binary Tree (Treap) + + --------------------------------------------------------------------- + Copyright (C) 2026 + Jonathan Fieldsend + Manuel Lopez-Ibanez + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + --------------------------------------------------------------------- + +*************************************************************************/ +#include +#include +#include +#include "sort.h" +#include "maxminclamp.h" + +#include "treap_archive.h" +#include "archiving_priv.h" + +#include "treap_iterator.h" + + +//--------------------------------------------------------------------------- + +bool +treap_archive_move(TreapArchive * restrict dst, TreapArchive * restrict src) +{ + memcpy(dst, src, sizeof(*src)); + memset(src, 0, sizeof(*src)); + free(src); + return true; +} + +void +treap_archive_init(TreapArchive *arch) +{ + treap_init(&arch->tree); + arch->total_size = 0; + arch->unique_size = 0; + arch->revision = 0; + arch->x_is_null = -1; // Don't know yet. +} + + +TreapArchive * +treap_archive_new(void) +{ + TreapArchive * arch = malloc(sizeof(*arch)); + if (arch != NULL) + treap_archive_init(arch); + + return arch; +} + +static inline TreapNode * +treap_archive_node_new(const double * restrict z, const void * restrict x) +{ + TreapNode * node = malloc(sizeof(*node)); + if (!node) + return NULL; + if (!void_list_init(&node->item, x)) { + free(node); + return NULL; + } + treap_node_init(node, z[0], z[1], node->item); + return node; +} + +bool +treap_archive_init_with_single_solution(TreapArchive * restrict arch, + const double * restrict z, const void * restrict x) +{ + // FIXME: Use a pool of nodes. + TreapNode *node = treap_archive_node_new(z, x); + if (node == NULL) + return false; + + assert(node->left == NULL); + assert(node->right == NULL); + treap_init_with_single_node(&arch->tree, node); + arch->unique_size = 1; + arch->total_size = node->item.n; + arch->revision = 1; + arch->x_is_null = (node->item.list == NULL); + return true; +} + +void +treap_archive_free_nodes(TreapArchive * arch) +{ + treap_free_nodes(&arch->tree); +} + +void +treap_archive_free(TreapArchive * arch) +{ + treap_archive_free_nodes(arch); + free(arch); +} + +size_t +treap_archive_total_size(const TreapArchive *arch) +{ + return arch->total_size; +} + +size_t +treap_archive_unique_size(const TreapArchive *arch) +{ + return arch->unique_size; +} + +bool +treap_archive_find_exact_vector(const TreapArchive * arch, const double * z) +{ + TreapNode *node = treap_find(&arch->tree, z[0]); + return (node && node->y == z[1]); +} + +VoidList * +treap_archive_find_solutions_of_vector(const TreapArchive * arch, const double * z) +{ + TreapNode *node = treap_find(&arch->tree, z[0]); + if (node && node->y == z[1]) + return &node->item; + return NULL; +} + + +/** + Returns true if z dominates every point in the tree, false otherwise or if z + exists in the tree or the tree is empty. + + Complexity: O(2 log n) +*/ +bool +treap_archive_dominated_by(const TreapArchive * arch, const double * z) +{ + const TreapNode * root = arch->tree.root; + if (arch->unique_size == 0) { + assert(root == NULL); + return false; + } + if (arch->unique_size == 1) { + return dominates_2d(z, (double [2]){ root->x, root->y }); + } + + const TreapNode * first = treap_minimum(root); + const TreapNode * last = treap_maximum(root); + assert(first != last); + // Even if they are both equal, first and last are not the same point. + return z[0] <= first->x && z[1] <= last->y; +} + +bool +treap_archive_dominates(const TreapArchive * arch, const double * z) +{ + return treap_dominates(&arch->tree, z); +} + +// FIXME: This visits every node, so it is O(log n). Is there a more efficient way? +static inline size_t +treap_archive_subtree_unique_size(const TreapNode *root, size_t * total_size) +{ + if (!root) + return 0; + *total_size += root->item.n; + return 1 + treap_archive_subtree_unique_size(root->left, total_size) + + treap_archive_subtree_unique_size(root->right, total_size); +} + +static inline void +treap_archive_update_size(TreapArchive *arch) +{ + arch->total_size = 0; + arch->unique_size = treap_archive_subtree_unique_size(arch->tree.root, &arch->total_size); +} + +/** + Insert new_node into self and return the nodes it displaces. + + Preconditions: + + - new_node is initialized and has a valid priority. + - new_node->left == NULL + - new_node->right == NULL + - new_node->z[0] is its key. + - its key does not already occur in self. + - no node in self weakly dominates new_node. + - self is a valid mutually-non-dominated minimization frontier. + + Postconditions: + + - new_node belongs to self. + - every node returned is dominated by new_node. + - no other node is removed from self. + - the returned nodes form a valid treap/frontier. +*/ +static inline TreapArchive * +treap_archive_insert_node_and_displace(TreapArchive *self, TreapNode *node) +{ + assert(node != NULL); + + TreapArchive * displaced = treap_archive_new(); + if (!displaced) + return NULL; + + displaced->x_is_null = self->x_is_null; + displaced->tree.root = treap_insert_and_displace(&self->tree, node); + treap_archive_update_size(displaced); // O(n) ! Is there a more efficient way? + self->unique_size++; + self->unique_size -= displaced->unique_size; + self->total_size += node->item.n; + self->total_size -= displaced->total_size; + self->revision++; + return displaced; +} + +// FIXME: Free the displaced nodes as soon as they are found. Do not create a +// temporary Treap. +static inline void +treap_archive_insert_node_and_discard(TreapArchive *self, TreapNode *node) +{ + TreapArchive * displaced = treap_archive_insert_node_and_displace(self, node); + treap_archive_free(displaced); +} + + +/** + Insert new solution (z,x) and get the displaced treap nodes if displaced != NULL. + + Returns one of the values of archive_insert_result_t. +*/ +// FIXME: Split into check nocheck variants. +int +treap_archive_add(TreapArchive * restrict self, const double * restrict z, + const void * restrict x, bool check, TreapArchive ** restrict displaced) +{ + // FIXME: How to add to an already built tree of displaced? + if (displaced) { + assert(*displaced == NULL); + *displaced = NULL; + } + + TreapNode * leq = treap_find_le(&self->tree, z[0]); + if (leq && all_equal_double((double [2]){leq->x, leq->y}, z, 2)) { + /* Exact duplicates cannot change the front geometry; append them to + the existing treap node and skip all dominance work. */ + int res = wrong_x_is_null(&self->x_is_null, x); + if (res) + return res; + self->total_size++; + if (!void_list_append(&leq->item, x)) + return ARCHIVE_INSERT_MEMORY_ERROR; + return ARCHIVE_INSERT_DUPLICATED; // Accepted + } + + if (check && leq && leq->y <= z[1]) + return ARCHIVE_INSERT_REJECTED; + + int res = wrong_x_is_null(&self->x_is_null, x); + if (res) + return res; + + TreapNode *new_node = treap_archive_node_new(z, x); + if (new_node == NULL) + return ARCHIVE_INSERT_MEMORY_ERROR; // ERROR + + if (displaced) { + *displaced = treap_archive_insert_node_and_displace(self, new_node); + if (*displaced == NULL) + return ARCHIVE_INSERT_MEMORY_ERROR; + } else + treap_archive_insert_node_and_discard(self, new_node); + + return ARCHIVE_INSERT_ACCEPTED; // Accepted +} + +int +treap_archive_add_discard(TreapArchive * restrict self, + const double * restrict z, const void * restrict x) +{ + return treap_archive_add(self, z, x, /*check=*/true, /*displaced=*/NULL); +} + +/** + Insert the elements of incoming into arch and return displaced nodes in incoming. + It does not allocate any new memory! + + This function assumes that: + + - no vector in arch can weakly dominated a vector in incoming. + - no vector in incoming appears in arch. +*/ +int +treap_archive_insert_displaced(TreapArchive *arch, TreapArchive *incoming, bool *all_displaced) +{ + if (arch->x_is_null == -1) + arch->x_is_null = incoming->x_is_null; + // Failing this assert is a programming error and should never happen. + assert(arch->x_is_null == incoming->x_is_null); + + arch->unique_size += incoming->unique_size; + arch->total_size += incoming->total_size; + TreapNode * displaced; + // FIXME: Regenerate all priorities of incoming nodes using arch->tree so + // that we do not need to pass arch->tree to the union. + arch->tree.root = treap_frontier_union(&arch->tree, incoming->tree.root, arch->tree.root, &displaced, all_displaced); + // Reuse incoming for displaced. + incoming->tree.root = displaced; + treap_archive_update_size(incoming); // O(n) ! Is there a more efficient way? + arch->unique_size -= incoming->unique_size; + arch->total_size -= incoming->total_size; + arch->revision++; + return ARCHIVE_INSERT_ACCEPTED; +} + +bool +treap_archive_has_x_values(const TreapArchive *arch) +{ + // -1 or 1 means that there are not x values to retrieve. + return arch->x_is_null == 0; +} + +static inline void +treap_archive_node_get_x_values(const TreapNode *node, const void *** x_list) +{ + if (!node) + return; + + treap_archive_node_get_x_values(node->left, x_list); + void_list_get_x_values(&node->item, x_list); + treap_archive_node_get_x_values(node->right, x_list); +} + +void +treap_archive_get_x_values(const TreapArchive *arch, const void ** x_list) +{ + if (!treap_archive_has_x_values(arch)) + return; + treap_archive_node_get_x_values(arch->tree.root, &x_list); +} + +static inline void +treap_archive_node_get_contents(const TreapNode *node, double ** z_list, const void *** x_list) +{ + if (!node) + return; + + treap_archive_node_get_contents(node->left, z_list, x_list); + // Duplicate z for each x + for (size_t i = 0; i < node->item.n; i++, *z_list += 2) + memcpy(*z_list, (double [2]){ node->x, node->y}, 2 * sizeof **z_list); + + if (x_list != NULL) + void_list_get_x_values(&node->item, x_list); + treap_archive_node_get_contents(node->right, z_list, x_list); +} + +/** + z_list must be an array of double allocated with sufficient space. Vectors + are copied to this array. + + x_list must be of (void *). Only the pointers are copied. +*/ +void +treap_archive_get_contents(const TreapArchive *arch, double * z_list, const void ** x_list) +{ + treap_archive_node_get_contents(arch->tree.root, &z_list, + treap_archive_has_x_values(arch) ? &x_list : NULL); +} + +/** + z_list must be an array of pointers to 'const double' of length treap_unique_size(t). + + That is, z_list does not contain the actual data but just pointers into the + Treap. The data is not owned by the Treap so it may get deleted if not used + or copied immediately. +*/ +void +treap_archive_get_unique_vectors(const TreapArchive *arch, const double ** z_list) +{ + treap_get_unique_vectors(&arch->tree, z_list); +} diff --git a/c/treap_archive.h b/c/treap_archive.h new file mode 100644 index 000000000..5cb04ff6d --- /dev/null +++ b/c/treap_archive.h @@ -0,0 +1,74 @@ +/************************************************************************* + + treap-archive.h: 2D Archive based on Randomised Binary Tree (Treap) + + --------------------------------------------------------------------- + Copyright (C) 2026 + Jonathan Fieldsend + Manuel Lopez-Ibanez + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + --------------------------------------------------------------------- + +*************************************************************************/ +#ifndef TREAP_ARCHIVE_H +#define TREAP_ARCHIVE_H + +#include "archiving.h" + +typedef struct TreapArchive TreapArchive; +typedef struct TreapIterator TreapIterator; + +// FIXME: We have to include this here because archiver.h needs it. +#include "solutions_list.h" + +/** + The treap stores the 2D vectors. Optional x-values are stored in this VoidList. +*/ +typedef VoidList TreapItem; + +static inline void +treap_item_free(TreapItem item) +{ + void_list_free(&item); +} + +#include "treap2d.h" + +typedef struct TreapArchive { + Treap tree; + size_t unique_size; // Number of nodes in the tree, i.e., unique z vectors. + size_t total_size; // Number of x values (including duplicates of z vectors). + signed char x_is_null; + uint64_t revision; // To detect modifications of the tree while iterating. +} TreapArchive; + + +TreapArchive * treap_archive_new(void); +void treap_archive_init(TreapArchive *arch); +int treap_archive_add(TreapArchive * restrict self, const double * restrict z, + const void * restrict x, bool check, TreapArchive ** restrict displaced); +int treap_archive_add_discard(TreapArchive * restrict self, const double * restrict z, const void * restrict x); +int treap_archive_insert_displaced(TreapArchive *arch, TreapArchive *incoming, bool *all_displaced); +size_t treap_archive_total_size(const TreapArchive *arch); +size_t treap_archive_unique_size(const TreapArchive *arch); +void treap_archive_get_contents(const TreapArchive *arch, double * z_list, const void ** x_list); +void treap_archive_get_unique_vectors(const TreapArchive *arch, const double ** z_list); +void treap_archive_get_x_values(const TreapArchive *arch, const void ** x_list); +bool treap_archive_has_x_values(const TreapArchive *arch); +void treap_archive_free_nodes(TreapArchive * arch); +void treap_archive_free(TreapArchive * arch); +bool treap_archive_dominated_by(const TreapArchive * arch, const double * z); +bool treap_archive_dominates(const TreapArchive * arch, const double * z); +bool treap_archive_find_exact_vector(const TreapArchive * arch, const double * z); +VoidList * treap_archive_find_solutions_of_vector(const TreapArchive * arch, const double * z); + +typedef struct TreapIterator TreapIterator; +void treap_archive_iter_free(TreapIterator *it); +TreapIterator * treap_archive_iter_new(const TreapArchive *arch); +int treap_archive_iter_next(TreapIterator * restrict it, const double ** restrict z, const void ** restrict x); + +#endif // TREAP_ARCHIVE_H diff --git a/c/treap_iterator.h b/c/treap_iterator.h new file mode 100644 index 000000000..f70f0d97e --- /dev/null +++ b/c/treap_iterator.h @@ -0,0 +1,126 @@ +#ifndef TREAP_ITERATOR_H +#define TREAP_ITERATOR_H + +#define TREAP_ITERATOR_STACK_MAX 64 + +typedef struct TreapIterator { + const TreapNode * stack[TREAP_ITERATOR_STACK_MAX]; + uint32_t stack_size; + const TreapNode * current_node; + uint32_t solution_index; + // To detect modifications of the tree while iterating. + uint64_t revision; + const TreapArchive *arch; +} TreapIterator; + + +static bool +treap_archive_iter_push(TreapIterator * restrict it, const TreapNode * restrict node) +{ + assert(it->stack_size < TREAP_ITERATOR_STACK_MAX); + if (it->stack_size == TREAP_ITERATOR_STACK_MAX) + return false; + + it->stack[it->stack_size++] = node; + return true; +} + +static const TreapNode * +treap_archive_iter_pop(TreapIterator *it) +{ + return it->stack_size ? it->stack[--it->stack_size] : NULL; +} + + +void +treap_archive_iter_free(TreapIterator *it) +{ + if (it != NULL) + free(it); +} + +void +treap_archive_iter_init(TreapIterator * restrict it, const TreapArchive * restrict arch) +{ + it->stack_size = 0; + it->current_node = NULL; + it->solution_index = 0; + it->revision = arch->revision; + it->arch = arch; + /* Push the entire left spine of the tree. The smallest key will then be on + top of the stack. */ + const TreapNode *node = arch->tree.root; + while (node != NULL) { + if (!treap_archive_iter_push(it, node)) { + // FIXME: Signal an error! + return; + } + node = node->left; + } +} + +TreapIterator * +treap_archive_iter_new(const TreapArchive *arch) +{ + if (arch == NULL) + return NULL; + + TreapIterator *it = malloc(sizeof(*it)); + if (it == NULL) + return NULL; + + treap_archive_iter_init(it, arch); + return it; +} + +/** + 1 next element returned + 0 end of iteration + -1 error +*/ +int +treap_archive_iter_next(TreapIterator * restrict it, const double ** restrict z, const void ** restrict x) +{ + if (it == NULL || z == NULL || it->revision != it->arch->revision) + return -1; + + // Continue emitting solutions from the node most recently popped. + if (it->current_node != NULL) { + const VoidList *item = &it->current_node->item; + if (it->solution_index < item->n) { + *z = &it->current_node->x; + if (x) + *x = item->list[it->solution_index]; + it->solution_index++; + return 1; + } + } + + // Find the next node in in-order traversal. + const VoidList *item; + do { + if (it->stack_size == 0) { + *z = NULL; + if (x) + *x = NULL; + return 0; + } + + it->current_node = treap_archive_iter_pop(it); + const TreapNode *node = it->current_node->right; + while (node != NULL) { + if (!treap_archive_iter_push(it, node)) + return -1; + node = node->left; + } + item = &it->current_node->item; + } while (item->n == 0); + + *z = &it->current_node->x; + if (x) + *x = item->list[0]; + it->solution_index = 1; + return 1; +} + +#endif // TREAP_ITERATOR_H diff --git a/python/Makefile b/python/Makefile index 322850c70..1ffadd103 100644 --- a/python/Makefile +++ b/python/Makefile @@ -34,6 +34,9 @@ doc: fastdoc: $(MAKE) -C doc html-noplot +fulldoc: + $(MAKE) -C doc clean html + clean: $(MAKE) -C doc clean $(MAKE) -C src/moocore/libmoocore/ clean diff --git a/python/benchmarks/bench_archive_unbounded.py b/python/benchmarks/bench_archive_unbounded.py new file mode 100644 index 000000000..fd6d27506 --- /dev/null +++ b/python/benchmarks/bench_archive_unbounded.py @@ -0,0 +1,124 @@ +from bench import Bench + +import itertools +import numpy as np +import matplotlib.pyplot as plt + +import moocore + +from moarchiving import get_mo_archive as moarch_get_mo_archive + +# FIXME: How to test both float and Fractions? +moarch_get_mo_archive.hypervolume_computation_float_type = float +moarch_get_mo_archive.hypervolume_final_float_type = float + +datasets = { + "hypercube-uniform": lambda max_size, n_obj, rng: moocore.generate_sequence( + max_size, n_obj, method="cube", seed=rng + ), + "hypercube-uniform-rep4": lambda max_size, n_obj, rng: ( + moocore.generate_sequence( + int(max_size / 4), n_obj, method="cube", seed=rng, n_rep=max_size + ) + ), + "hypercube-ndsort=1": lambda max_size, n_obj, rng: ( + moocore.generate_sequence( + max_size, n_obj, method="cube", ndsort=1, seed=rng + ) + ), + "hypercube-ndsort=-1": lambda max_size, n_obj, rng: ( + moocore.generate_sequence( + max_size, n_obj, method="cube", ndsort=-1, seed=rng + ) + ), + "within-concave-sphere": lambda max_size, n_obj, rng: ( + moocore.generate_sequence( + max_size, n_obj, method="sphere", r_min=0.5, seed=rng + ) + ), + "inverted-simplex": lambda max_size, n_obj, rng: moocore.generate_ndset( + max_size, n_obj, method="inverted-simplex", seed=rng + ), + "glas2017,c=0.5": lambda max_size, n_obj, rng: moocore.generate_sequence( + max_size, n_obj, method="glas2017", c_value=0.5, seed=rng + ), + "glas2017,c=2.0": lambda max_size, n_obj, rng: moocore.generate_sequence( + max_size, n_obj, method="glas2017", c_value=2.0, seed=rng + ), + "cliff-convex": lambda max_size, n_obj, rng: moocore.generate_ndset( + max_size, n_obj, method="cliff-convex", seed=rng + ), + "cliff-concave": lambda max_size, n_obj, rng: moocore.generate_ndset( + max_size, n_obj, method="cliff-concave", seed=rng + ), +} + +n_repeats = 10 + +title = "Adding a point to the archive" +file_prefix = "unbarchive_add" +max_exp = {2: 19, 3: 17, 4: 13} +print(f"Running benchmark: {title}") +names = datasets.keys() +n_objs = [2, 3, 4] +for n_obj, name in itertools.product(n_objs, names): + if n_obj == 2 and name.startswith("cliff"): + continue + test_archive_sizes = 2 ** np.arange(1, max_exp[n_obj]) + max_size = test_archive_sizes[-1] + + # This is a bit wasteful because we rebuild the archive from scratch for + # each size. On the other hand, it is probably less sensitive to noise. + setup = { + "moocore (UnboundedArchive)": lambda z: dict( + moa=moocore.UnboundedArchive(dim=n_obj), z=z + ), + "moarchiving (float)": lambda z: dict( + moa=moarch_get_mo_archive(n_obj=n_obj), z=z.tolist() + ), + } + benchmarks = { + "moocore (UnboundedArchive)": lambda moa, z: [moa.add(p) for p in z], + "moocore (filter_dominated)": lambda z: moocore.filter_dominated(z), + "moarchiving (float)": lambda moa, z: [moa.add(p) for p in z], + } + # Divide time by number of points, so that we report the average time per + # point for different input sizes. + normalise = ["moocore (UnboundedArchive)", "moarchiving (float)"] + + if n_obj == 2: + setup["moocore (NDTree)"] = lambda z: dict( + moa=moocore._archive.NDTreeArchive(dim=n_obj), z=z + ) + benchmarks["moocore (NDTree)"] = lambda moa, z: [moa.add(p) for p in z] + normalise.append("moocore (NDTree)") + + bench = Bench( + name=f"{name}-{n_obj}d", + n=test_archive_sizes, + bench=benchmarks, + setup=setup, + reps=1, + max_time=5, + ) + + for r in range(n_repeats): + # We generate the dataset outside setup because we want to do it once + # per repetition, not once per benchmark. + rng = np.random.default_rng(42 + r) + z = datasets[name](max_size, n_obj, rng) + bench(lambda n: z[:n]) + + bench.normalise_times(normalise) + + bench.plots( + file_prefix=file_prefix, + title=title, + log="xy", + xlabel="Number of points", + logx_base=2, + show_ci=True, + ) + +if "__file__" not in globals(): # Running interactively. + plt.show() diff --git a/python/doc/source/REFERENCES.bib b/python/doc/source/REFERENCES.bib index b752aab29..871755195 100644 --- a/python/doc/source/REFERENCES.bib +++ b/python/doc/source/REFERENCES.bib @@ -298,6 +298,42 @@ @article{HerWer1987tabucol pages = {345--351} } +@article{JasLus2018ndtree, + author = { Andrzej Jaszkiewicz and Thibaut Lust }, + title = {{ND}-Tree-Based Update: A Fast Algorithm for the Dynamic + Nondominance Problem}, + journal = {IEEE Transactions on Evolutionary Computation}, + year = 2018, + volume = 22, + number = 5, + pages = {778--791}, + doi = {10.1109/TEVC.2018.2799684}, + publisher = {IEEE}, + abstract = {In this paper, we propose a new method called ND-Tree-based + update (ND-Tree) for the dynamic nondominance problem, i.e., + the problem of online update of a Pareto archive composed of + mutually nondominated points. It uses a new ND-Tree data + structure in which each node represents a subset of points + contained in a hyperrectangle defined by its local + approximate ideal and nadir points. By building subsets + containing points located close in the objective space and + using basic properties of the local ideal and nadir points we + can efficiently avoid searching many branches in the + tree. ND-Tree may be used in multiobjective evolutionary + algorithms and other multiobjective metaheuristics to update + an archive of potentially nondominated points. We prove that + the proposed algorithm has sublinear time complexity under + mild assumptions. We experimentally compare ND-Tree to the + simple list, Quad-tree, and M-Front methods using artificial + and realistic benchmarks with up to ten objectives and show + that with this new method substantial reduction of the number + of point comparisons and computational time can be + obtained. Furthermore, we apply the method to the + nondominated sorting problem showing that it is highly + competitive to some recently proposed algorithms dedicated to + this problem.} +} + @article{Jen03, author = {M. T. Jensen}, title = {Reducing the run-time complexity of multiobjective {EA}s: The @@ -455,6 +491,16 @@ @article{ZitThiLauFon2003:tec doi = {10.1109/TEVC.2003.810758} } +@inproceedings{AraSei1989rst, + author = {Aragon, C. R. and Seidel, R. G.}, + title = {Randomized search trees}, + booktitle = {30th Annual Symposium on Foundations of Computer Science}, + year = 1989, + pages = {540--545}, + keywords = {Binary search trees, Treap, Dynamic programming}, + doi = {10.1109/SFCS.1989.63531} +} + @incollection{AugBadBroZit2009gecco, booktitle = {Proceedings of the Genetic and Evolutionary Computation Conference, GECCO 2009}, address = { New York, NY}, diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-3d-time.png new file mode 100644 index 000000000..28f646b70 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-4d-time.png new file mode 100644 index 000000000..1064a64eb Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-concave-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-3d-time.png new file mode 100644 index 000000000..2ab0549ee Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-4d-time.png new file mode 100644 index 000000000..050e83770 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-cliff-convex-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-2d-time.png new file mode 100644 index 000000000..5e3bd77de Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-3d-time.png new file mode 100644 index 000000000..535836971 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-4d-time.png new file mode 100644 index 000000000..6be53d7dd Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=0.5-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-2d-time.png new file mode 100644 index 000000000..2f0eb5f8d Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-3d-time.png new file mode 100644 index 000000000..d085d4bee Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-4d-time.png new file mode 100644 index 000000000..370d1d657 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-glas2017,c=2.0-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-2d-time.png new file mode 100644 index 000000000..0b02d336a Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-3d-time.png new file mode 100644 index 000000000..f64c207a9 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-4d-time.png new file mode 100644 index 000000000..73f7d3de6 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=-1-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-2d-time.png new file mode 100644 index 000000000..699652d92 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-3d-time.png new file mode 100644 index 000000000..bdf8b4e47 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-4d-time.png new file mode 100644 index 000000000..394c6c03d Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-ndsort=1-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-2d-time.png new file mode 100644 index 000000000..dda3b1fc4 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-3d-time.png new file mode 100644 index 000000000..85ba950d8 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-4d-time.png new file mode 100644 index 000000000..ae594ed2a Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-2d-time.png new file mode 100644 index 000000000..31d9cb50c Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-3d-time.png new file mode 100644 index 000000000..7e312007f Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-4d-time.png new file mode 100644 index 000000000..9e5ef2111 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-hypercube-uniform-rep4-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-2d-time.png new file mode 100644 index 000000000..e9967f981 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-3d-time.png new file mode 100644 index 000000000..a6e0de8b0 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-4d-time.png new file mode 100644 index 000000000..0a11b3276 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-inverted-simplex-4d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-2d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-2d-time.png new file mode 100644 index 000000000..e11342c3f Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-2d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-3d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-3d-time.png new file mode 100644 index 000000000..b0644dd24 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-3d-time.png differ diff --git a/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-4d-time.png b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-4d-time.png new file mode 100644 index 000000000..6565fddb0 Binary files /dev/null and b/python/doc/source/_static/bench/unbarchive_add_bench-within-concave-sphere-4d-time.png differ diff --git a/python/doc/source/_templates/autosummary/class.rst b/python/doc/source/_templates/autosummary/class.rst index 4e2358f8b..f7bda6399 100644 --- a/python/doc/source/_templates/autosummary/class.rst +++ b/python/doc/source/_templates/autosummary/class.rst @@ -10,7 +10,7 @@ .. autoclass:: {{ objname }} :members: :inherited-members: - :special-members: __call__ + :special-members: __call__, __len__, __iter__ .. minigallery:: {{ module }}.{{ objname }} {% for meth in methods %}{{ module }}.{{ objname }}.{{ meth }} {% endfor %} :add-heading: Gallery examples diff --git a/python/doc/source/benchmarks.inc.rst b/python/doc/source/benchmarks.inc.rst index 951d692b7..072c7d39b 100644 --- a/python/doc/source/benchmarks.inc.rst +++ b/python/doc/source/benchmarks.inc.rst @@ -260,6 +260,11 @@ The following plots benchmark the speed of computing the :ref:`epsilon indicator :width: 48% +Archiving Benchmarks +-------------------- + +See :ref:`Archiving Benchmarks `. + .. _BoFire: https://experimental-design.github.io/bofire/ .. _BoTorch: https://botorch.org/ diff --git a/python/doc/source/benchmarks_archiving.rst b/python/doc/source/benchmarks_archiving.rst new file mode 100644 index 000000000..14e83af02 --- /dev/null +++ b/python/doc/source/benchmarks_archiving.rst @@ -0,0 +1,162 @@ +:orphan: + +.. _benchmarks_archiving: + +#################### +Archiving Benchmarks +#################### + +The following plots compare the performance of archiving algorithms in `moocore`_ and `moarchiving`_. + +The sequences are generated using :func:`~moocore.generate_ndset` and :func:`~moocore.generate_sequence`. + +The source code for the benchmarks below can be found at https://github.com/multi-objective/moocore/tree/main/python/benchmarks . + +Unbounded Archiving +------------------- + +This benchmark compares :func:`moocore.UnboundedArchive` and +|moarchiving_get_mo_archive|_. Each archive is given the number of points in +the x-axis, one at a time. We measure the CPU time required to archive the +whole sequence of points and divide it by the number of points, which gives an +estimation of the time to process each point for different input sizes. In +addition, the plots show the total time required by +:func:`~moocore.filter_dominated` to process the input sequence. The time +required by :func:`~moocore.filter_dominated` can be seen as an upper-bound. + +2D +~~ + +|unbarchive_add_bench-glas2017,c=0.5-2d-time| |unbarchive_add_bench-glas2017,c=2.0-2d-time| + +|unbarchive_add_bench-hypercube-uniform-2d-time| |unbarchive_add_bench-hypercube-uniform-rep4-2d-time| + +|unbarchive_add_bench-hypercube-ndsort=-1-2d-time| |unbarchive_add_bench-hypercube-ndsort=1-2d-time| + +|unbarchive_add_bench-inverted-simplex-2d-time| |unbarchive_add_bench-within-concave-sphere-2d-time| + +.. |unbarchive_add_bench-glas2017,c=0.5-2d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=0.5-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-glas2017,c=2.0-2d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=2.0-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-inverted-simplex-2d-time| image:: _static/bench/unbarchive_add_bench-inverted-simplex-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-2d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-rep4-2d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-rep4-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=-1-2d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=-1-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=1-2d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=1-2d-time.png + :width: 49% + +.. |unbarchive_add_bench-within-concave-sphere-2d-time| image:: _static/bench/unbarchive_add_bench-within-concave-sphere-2d-time.png + :width: 49% + + + +3D +~~ + + +|unbarchive_add_bench-glas2017,c=0.5-3d-time| |unbarchive_add_bench-glas2017,c=2.0-3d-time| + +|unbarchive_add_bench-hypercube-uniform-3d-time| |unbarchive_add_bench-hypercube-uniform-rep4-3d-time| + +|unbarchive_add_bench-hypercube-ndsort=-1-3d-time| |unbarchive_add_bench-hypercube-ndsort=1-3d-time| + +|unbarchive_add_bench-inverted-simplex-3d-time| |unbarchive_add_bench-within-concave-sphere-3d-time| + +|unbarchive_add_bench-cliff-concave-3d-time| |unbarchive_add_bench-cliff-convex-3d-time| + +.. |unbarchive_add_bench-glas2017,c=0.5-3d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=0.5-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-glas2017,c=2.0-3d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=2.0-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-inverted-simplex-3d-time| image:: _static/bench/unbarchive_add_bench-inverted-simplex-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-3d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-rep4-3d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-rep4-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=-1-3d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=-1-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=1-3d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=1-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-within-concave-sphere-3d-time| image:: _static/bench/unbarchive_add_bench-within-concave-sphere-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-cliff-concave-3d-time| image:: _static/bench/unbarchive_add_bench-cliff-concave-3d-time.png + :width: 49% + +.. |unbarchive_add_bench-cliff-convex-3d-time| image:: _static/bench/unbarchive_add_bench-cliff-convex-3d-time.png + :width: 49% + +4D +~~ + +|unbarchive_add_bench-glas2017,c=0.5-4d-time| |unbarchive_add_bench-glas2017,c=2.0-4d-time| + +|unbarchive_add_bench-hypercube-uniform-4d-time| |unbarchive_add_bench-hypercube-uniform-rep4-4d-time| + +|unbarchive_add_bench-hypercube-ndsort=-1-4d-time| |unbarchive_add_bench-hypercube-ndsort=1-4d-time| + +|unbarchive_add_bench-inverted-simplex-4d-time| |unbarchive_add_bench-within-concave-sphere-4d-time| + +|unbarchive_add_bench-cliff-concave-4d-time| |unbarchive_add_bench-cliff-convex-4d-time| + +.. |unbarchive_add_bench-glas2017,c=0.5-4d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=0.5-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-glas2017,c=2.0-4d-time| image:: _static/bench/unbarchive_add_bench-glas2017,c=2.0-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-inverted-simplex-4d-time| image:: _static/bench/unbarchive_add_bench-inverted-simplex-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-4d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-uniform-rep4-4d-time| image:: _static/bench/unbarchive_add_bench-hypercube-uniform-rep4-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=-1-4d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=-1-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-hypercube-ndsort=1-4d-time| image:: _static/bench/unbarchive_add_bench-hypercube-ndsort=1-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-within-concave-sphere-4d-time| image:: _static/bench/unbarchive_add_bench-within-concave-sphere-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-cliff-concave-4d-time| image:: _static/bench/unbarchive_add_bench-cliff-concave-4d-time.png + :width: 49% + +.. |unbarchive_add_bench-cliff-convex-4d-time| image:: _static/bench/unbarchive_add_bench-cliff-convex-4d-time.png + :width: 49% + + + +.. _moarchiving: https://cma-es.github.io/moarchiving/ +.. _moocore: https://multi-objective.github.io/moocore/python/ + +.. The :func: role expects a Python cross-reference, not an arbitrary URL. + Since reStructuredText does not support nested inline markup, use a + substitution to combine code formatting with an external hyperlink. + See https://www.docutils.org/FAQ.html#is-nested-inline-markup-possible + +.. |moarchiving_get_mo_archive| replace:: :code:`moarchiving.get_mo_archive` +.. _moarchiving_get_mo_archive: https://cma-es.github.io/moarchiving/moarchiving-apidocs/moarchiving.get_archive.html#get_mo_archive diff --git a/python/doc/source/conf.py b/python/doc/source/conf.py index 4a80c3dd6..bb2c69c6d 100644 --- a/python/doc/source/conf.py +++ b/python/doc/source/conf.py @@ -97,6 +97,7 @@ # If true, '()' will be appended to :func: etc. cross-reference text. add_function_parentheses = True +autodoc_inherit_docstrings = True autodoc_typehints = "none" # Conflicts with sphinx_autodoc_typehints typehints_document_rtype = True typehints_use_rtype = False @@ -108,6 +109,8 @@ "ArrayLike": ":py:data:`~numpy.typing.ArrayLike`", "NDArray": ":py:data:`~numpy.typing.NDArray`", } +# Hide the class/module prefixes for objects in TOC entries (sidebar) +toc_object_entries_show_parents = "hide" nitpicky = True nitpick_ignore = { @@ -125,7 +128,6 @@ bibtex_bibliography_id = "bibliography-{bibliography_count}" bibtex_footbibliography_id = "footbibliography-{footbibliography_count}" - # FIXME: This is not working! logger = sphinx.util.logging.getLogger("sphinx.ext.autodoc") members_to_watch = ["class", "exception", "function", "method", "attribute"] diff --git a/python/doc/source/reference/functions.archiving.rst b/python/doc/source/reference/functions.archiving.rst new file mode 100644 index 000000000..5dd46c607 --- /dev/null +++ b/python/doc/source/reference/functions.archiving.rst @@ -0,0 +1,18 @@ +.. _archiving: + +========================== +Archiving Pareto frontiers +========================== + +.. currentmodule:: moocore + +.. _unbounded_archiving: + +Unbounded Archiving of a Nondominated Set +========================================= + +.. autosummary:: + :toctree: generated/ + + UnboundedArchive + BaseUnboundedArchive diff --git a/python/doc/source/reference/index.rst b/python/doc/source/reference/index.rst index c28449611..943651b1e 100644 --- a/python/doc/source/reference/index.rst +++ b/python/doc/source/reference/index.rst @@ -29,3 +29,4 @@ indentation. functions.dominance functions.metrics functions.eaf + functions.archiving diff --git a/python/doc/source/whatsnew/index.rst b/python/doc/source/whatsnew/index.rst index ea58ce81e..3844e36ed 100644 --- a/python/doc/source/whatsnew/index.rst +++ b/python/doc/source/whatsnew/index.rst @@ -12,6 +12,8 @@ Version 0.4.0 - :func:`~moocore.is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation. - New shapes ``"concave-cliff"`` and ``"convex-cliff"`` added to :func:`~moocore.generate_ndset`. - New function :func:`~moocore.generate_sequence` to generate sequences of dominated and nondominated points. +- New unbounded nondominated archive :func:`~moocore.UnboundedArchive`. +- New archiving benchmarks :ref:`benchmarks_archiving`. Version 0.3.2 (11/07/2026) diff --git a/python/src/moocore/__init__.py b/python/src/moocore/__init__.py index d07f559f0..3808217f9 100644 --- a/python/src/moocore/__init__.py +++ b/python/src/moocore/__init__.py @@ -37,6 +37,10 @@ ) from ._generate import generate_ndset, generate_sequence +from ._archive import UnboundedArchive + +# The redundant alias silences a ruff warning. +from ._archive import BaseUnboundedArchive as BaseUnboundedArchive from importlib.metadata import version as _metadata_version @@ -49,6 +53,7 @@ "Hypervolume", "ReadDatasetsError", "RelativeHypervolume", + "UnboundedArchive", "any_dominated", "apply_within_sets", "avg_hausdorff_dist", diff --git a/python/src/moocore/_archive.py b/python/src/moocore/_archive.py new file mode 100644 index 000000000..423d06fee --- /dev/null +++ b/python/src/moocore/_archive.py @@ -0,0 +1,696 @@ +from abc import ABC, abstractmethod +from collections.abc import Callable, Sequence, Iterator +from numpy.typing import ArrayLike # For type hints +from typing import Any, ClassVar + +import numpy as np + +from ._docsubstitute import DocSubstitute +from ._utils import _as_1d_float_array, _as_2d_float_array, _parse_maximise + +## The CFFI library is used to create C bindings. +from ._libmoocore import lib, ffi # ty: ignore[unresolved-import] + + +def _remember_handle(handles, what): + handles[what] = what + + +class BaseArchive(ABC): + _lib_archive_free: ClassVar[Callable[..., Any]] + + def __init__(self, dim: int, maximise: bool | Sequence[bool]): + if dim < 2 or dim > lib.MOOCORE_DIMENSION_MAX: + raise ValueError( + f"'dim' cannot be smaller than 2 or larger than {lib.MOOCORE_DIMENSION_MAX}" + ) + + self._dim = dim + self._c_archive = self._setup_c_archive() + if self._c_archive == ffi.NULL: + raise MemoryError + # ffi.gc() ensures that the pointer is registered with the GC and the C + # destructor called. + self._c_archive = ffi.gc(self._c_archive, self._lib_archive_free) + + self.maximise = np.asarray( + _parse_maximise(maximise, self._dim), dtype=bool + ) + # -1: maximize, 1: minimize + self._minmax = np.where(self.maximise, -1.0, 1.0) + + # Python objects stored in C will be freed by Python unless we keep a + # reference to them in Python. This set keeps those objects alive until + # we decide to free or move them to another front. + self._x_handles = {} + # FIXME: Using a set() instead of a dictionary is probably faster. We + # can store the output of ffi.new_handle() directly in the set. + # However, when retrieving handles from C, two handles may compare + # equal, but the one in self._x_handles is the one that should be kept + # alive. + self.displaced_p = ffi.NULL + + @property + def dim(self) -> int: + """Return the dimension of the points (objective vectors) stored in the archive.""" + return self._dim + + @abstractmethod + def _setup_c_archive(self): + pass + + def _objectives_to_c(self, z: ArrayLike): + z = _as_1d_float_array(z, dim=self.dim, name="z") + # We do not want to modify the original array, so we create a copy here. + z = z * self._minmax + return ffi.from_buffer("double []", z) + + def _c_to_objectives(self, z) -> np.ndarray: + z = np.frombuffer(ffi.buffer(z, ffi.sizeof("double") * self.dim)) + z = z * self._minmax + z.setflags(write=False) + return z + + +class BaseUnboundedArchive(BaseArchive): + """Interface implemented by objects returned by :func:`UnboundedArchive`. + + This class defines the common interface for ``BaseUnboundedArchive`` + objects. It is not intended to be instantiated directly. See Examples in + :func:`UnboundedArchive`. + + """ + + # Derived classes must define these: + _lib_archive_add_and_get_displaced: ClassVar[Callable[..., Any]] + _lib_archive_dominated_by: ClassVar[Callable[..., Any]] + _lib_archive_dominates: ClassVar[Callable[..., Any]] + _lib_archive_find_exact_vector: ClassVar[Callable[..., Any]] + _lib_archive_find_x_values_of_vector: ClassVar[Callable[..., Any]] + _lib_archive_free: ClassVar[Callable[..., Any]] + _lib_archive_get_contents: ClassVar[Callable[..., Any]] + _lib_archive_get_unique_vectors: ClassVar[Callable[..., Any]] + _lib_archive_has_x_values: ClassVar[Callable[..., Any]] + _lib_archive_new: ClassVar[Callable[..., Any]] + _lib_archive_total_size: ClassVar[Callable[..., Any]] + _lib_archive_unique_size: ClassVar[Callable[..., Any]] + _lib_archive_iter_new: ClassVar[Callable[..., Any]] + _lib_archive_iter_free: ClassVar[Callable[..., Any]] + _lib_archive_iter_next: ClassVar[Callable[..., Any]] + _lib_displaced_free: ClassVar[Callable[..., Any]] + _lib_displaced_get_unique_vectors: ClassVar[Callable[..., Any]] + _lib_displaced_get_x_values: ClassVar[Callable[..., Any]] + _lib_displaced_total_len: ClassVar[Callable[..., Any]] + _lib_displaced_unique_len: ClassVar[Callable[..., Any]] + + def __init__(self, dim: int, maximise: bool | Sequence[bool]): + super().__init__(dim=dim, maximise=maximise) + self.all_displaced_p = ffi.new("bool *") + + def __contains__(self, z: ArrayLike) -> bool: + """Return true if the given objective vector is archived.""" + z = self._objectives_to_c(z) + return self._lib_archive_find_exact_vector(self._c_archive, z) + + def __iter__(self) -> Iterator[tuple[np.ndarray, Any]]: + """Return z vectors and, if available, their associated x object one at a time. + + For efficiency, iterating does not make copies of the vectors nor the + objects. However, the vectors are marked as read-only, thus + modifying them requires that the user makes a copy. + + Modifying the archive while iterating raises a :exc:`RuntimeError`. + """ + if len(self) == 0: + return + + it = self._lib_archive_iter_new(self._c_archive) + if it == ffi.NULL: + raise MemoryError("failed to create iterator") + it = ffi.gc(it, self._lib_archive_iter_free) + z = ffi.new("const double **") + rc = 1 + if self.has_x_values(): + x = ffi.new("const void **") + while True: + rc = self._lib_archive_iter_next(it, z, x) + if rc != 1: + break + yield self._c_to_objectives(z[0]), ffi.from_handle(x[0]) + else: + while True: + rc = self._lib_archive_iter_next(it, z, ffi.NULL) + if rc != 1: + break + yield self._c_to_objectives(z[0]), None + + if rc < 0: + raise RuntimeError( + "invalid 'BaseUnboundedArchive' iterator, did you modify the archive?" + ) + + def __len__(self) -> int: + """Return the total number of points stored in the archive (including duplicates).""" + return self._lib_archive_total_size(self._c_archive) + + def __repr__(self) -> str: + z, x = self.get_contents() + if x is None: + return str(z) + return f"{z}\nx = {x}" + + def unique_len(self) -> int: + """Return the unique number of points stored in the archive.""" + return self._lib_archive_unique_size(self._c_archive) + + def has_x_values(self) -> bool: + """Return true if the archive stores 'x' values.""" + return self._lib_archive_has_x_values(self._c_archive) + + def dominates(self, z: ArrayLike) -> bool: + """Return true if a point in the archive dominates :code:`z`.""" + z = self._objectives_to_c(z) + return self._lib_archive_dominates(self._c_archive, z) + + # FIXME: Should this be called is_dominated_by, so you write arch.is_dominated_by(z)? + def dominated_by(self, z: ArrayLike) -> bool: + """Return true if all points in the archive are dominated by :code:`z`.""" + z = self._objectives_to_c(z) + return self._lib_archive_dominated_by(self._c_archive, z) + + def add(self, z: ArrayLike, x: Any = None) -> bool: + """Add a point if it is not (weakly-)dominated by the archive. + + Points dominated by the added point are discarded. With + :code:`allow_duplicates=True`, weakly-dominated points are + accepted. Otherwise, they are rejected. + + Returns + ------- + ``True`` if added or ``False`` if rejected. + """ + accepted, _ = self._add_with_displaced_raw( + z, x, check=True, discard=True + ) + return accepted + + def add_many( + self, z: ArrayLike, x: Sequence[Any] | None = None + ) -> list[bool]: + """Add several points sequentially and return acceptance as a boolean list.""" + # FIXME: call is_nondominated(points) first. It is much faster than adding + # them to the archive one by one. + z = _as_2d_float_array(z, dim=self.dim) + if x is None: + return [self.add(point) for point in z] + + x = list(x) + if len(x) != len(z): + raise ValueError( + "'x' must be either None or a list with one value per input vector" + ) + return [self.add(point, x_value) for point, x_value in zip(z, x)] + + # FIXME: This should be called _add_and_get_displaced() + def _add_with_displaced_raw( + self, z: ArrayLike, x: Any, check: bool, discard: bool + ): + """Insert an entry and return members displaced into the next Pareto front.""" + # The C code copies z but not x. + z_c = self._objectives_to_c(z) + x_c = ffi.NULL if x is None else ffi.new_handle(x) + + # Make sure _lib_archive_add_and_get_displaced() starts from an empty + # list of displaced. + self.displaced_p[0] = ffi.NULL + # FIXME: Add option to discard removed points (so no displaced). + accepted = self._lib_archive_add_and_get_displaced( + self._c_archive, z_c, x_c, check, self.displaced_p + ) + if accepted == lib.ARCHIVE_INSERT_REJECTED: + return False, None + + if ( + accepted != lib.ARCHIVE_INSERT_ACCEPTED + and accepted != lib.ARCHIVE_INSERT_DUPLICATED + ): + if accepted == lib.ARCHIVE_INSERT_X_MUST_BE_NOT_NULL: + raise ValueError("'x' cannot be None") + elif accepted == lib.ARCHIVE_INSERT_X_MUST_BE_NULL: + raise ValueError("'x' must be None") + elif accepted == lib.ARCHIVE_INSERT_MEMORY_ERROR: + raise MemoryError + else: + raise NotImplementedError("unknown error") + + if x is not None: + _remember_handle(self._x_handles, x_c) + + # FIXME: if x is None and discard: return True, None + displaced = self.displaced_p[0] + if displaced == ffi.NULL: + return True, None + + displaced = ffi.gc(displaced, self._lib_displaced_free) + # FIXME: We know that displaced is non-zero, so the checks are wasted. + return True, self._process_displaced(displaced, discard=discard) + + def get_contents(self) -> tuple[np.ndarray, list[Any] | None]: + """Return the full contents of the archive. + + Returns + ------- + Tuple of objective vectors (as a Numpy array) and a list of x values (or ``None`` if ``not self.has_x_values()``) + + """ + return self._c_get_contents(self._c_archive, self.dim, self._minmax) + + def get_unique_vectors(self) -> np.ndarray: + """Return the unique objective vectors stored in the archive.""" + return self._c_get_unique_vectors( + self._c_archive, self.dim, self._minmax + ) + + def get_x_values(self) -> list[Any]: + """Return the x-values stored in the archive.""" + return self._c_get_x_values(self._c_archive) + + def get_x_values_of_vector(self, z: ArrayLike) -> list[Any]: + """Return the x-values of vector ``z``. + + Returns + ------- + A list of objects or an empty list if there are none or ``z`` is not in the archive. + + """ + z = self._objectives_to_c(z) + x_list = self._lib_archive_find_x_values_of_vector(self._c_archive, z) + if x_list == ffi.NULL: + return [] + return [ffi.from_handle(x_list.list[i]) for i in range(x_list.n)] + + def _get_x_handles_from_displaced(self, displaced): + if not self.has_x_values(): + return [] + + x_list = ffi.new("void *[]", self._lib_displaced_total_len(displaced)) + self._lib_displaced_get_x_values(displaced, x_list) + return x_list + + def _process_displaced(self, displaced, discard): + # FIXME: How can displaced be None? + # FIXME: How can displaced have len == 0? displaced should be NULL in that case. + if ( + displaced is None + or displaced == ffi.NULL + or self._lib_displaced_unique_len(displaced) == 0 + ): + return None + + handles = self._get_x_handles_from_displaced(displaced) + if discard: + for k in handles: + del self._x_handles[k] + return None + else: + x_handles = {k: self._x_handles.pop(k) for k in handles} + return displaced, x_handles + + def _add_displaced_raw(self, displaced): + if displaced is None: + return None + displaced, x_handles = displaced + displaced, all_displaced = self._lib_archive_insert_displaced(displaced) + if all_displaced: + assert self._lib_displaced_unique_len(displaced) > 0 + # Swap handles + x_handles, self._x_handles = self._x_handles, x_handles + return displaced, x_handles + + self._x_handles.update(x_handles) + return self._process_displaced(displaced, discard=False) + + @classmethod + def _c_get_contents(cls, c_tree, dim, minmax): + size = cls._lib_archive_total_size(c_tree) + z_list = ffi.new("double []", dim * size) + if cls._lib_archive_has_x_values(c_tree): + x_list = ffi.new("void *[]", size) + cls._lib_archive_get_contents(c_tree, z_list, x_list) + x_list = [ffi.from_handle(x) for x in x_list] + else: + cls._lib_archive_get_contents(c_tree, z_list, ffi.NULL) + x_list = [] + + z_list = np.frombuffer(ffi.buffer(z_list)).reshape(size, dim) + z_list *= minmax + return z_list, x_list + + @classmethod + def _c_get_unique_vectors(cls, c_archive, dim, minmax): + # FIXME: It would be faster to return a single block of memory, but we + # need a new function to get the pointer addresses. + n = cls._lib_archive_unique_size(c_archive) + if n == 0: + return [] + z_list = ffi.new("double *[]", n) + cls._lib_archive_get_unique_vectors(c_archive, z_list) + z_array = np.vstack( + [ + np.frombuffer(ffi.buffer(z, ffi.sizeof("double") * dim)) + for z in z_list + ] + ) + z_array *= minmax + return z_array + + @abstractmethod + def _lib_archive_insert_displaced(self, displaced) -> tuple[Any, bool]: + pass + + @staticmethod + @abstractmethod + def _c_get_x_values(c_archive) -> list[Any]: + pass + + +class NDTreeArchive(BaseUnboundedArchive): + # C functions + _lib_archive_add_and_get_displaced = lib.ndtree_add + _lib_archive_dominated_by = lib.ndtree_dominated_by + _lib_archive_dominates = lib.ndtree_dominates + _lib_archive_find_exact_vector = lib.ndtree_find_exact_vector + _lib_archive_find_x_values_of_vector = lib.ndtree_find_solutions_of_vector + _lib_archive_free = lib.ndtree_free + _lib_archive_get_contents = lib.ndtree_get_contents + _lib_archive_get_unique_vectors = lib.ndtree_get_unique_vectors + _lib_archive_has_x_values = lib.ndtree_has_x_values + _lib_archive_new = lib.ndtree_new + _lib_archive_total_size = lib.ndtree_total_size + _lib_archive_unique_size = lib.ndtree_unique_size + _lib_archive_iter_new = lib.ndtree_iter_new + _lib_archive_iter_free = lib.ndtree_iter_free + _lib_archive_iter_next = lib.ndtree_iter_next + _lib_displaced_free = lib.flex_bucket_free + _lib_displaced_get_unique_vectors = lib.flex_bucket_get_unique_vectors + _lib_displaced_get_x_values = lib.flex_bucket_get_x_values + _lib_displaced_total_len = lib.flex_bucket_total_len + _lib_displaced_unique_len = lib.flex_bucket_unique_len + + def __init__( + self, + dim: int, + maximise: bool | Sequence[bool] = False, + max_children: int = 0, + max_bucket_size: int = 0, + allow_duplicates: bool = True, + ): + + self.max_children, self.max_bucket_size, self.allow_duplicates = ( + self._get_default_config( + dim, max_children, max_bucket_size, allow_duplicates + ) + ) + super().__init__(dim=dim, maximise=maximise) + # We cache these pointers here so we don't keep allocating and freeing memory. + self.displaced_p = ffi.new("FlexBucket **") + + # This is a @classmethod so we can use it from NDSForestNDTree + @classmethod + def _get_default_config( + cls, + dim, + max_children: int, + max_bucket_size: int, + allow_duplicates: bool, + ): + # FIXME: This duplicates the logic in C, it would be better to call a C + # function. + max_children = max(dim + 1, max_children) + if max_children >= 255: + raise ValueError("'max_children' must be smaller than 255") + + if max_bucket_size < max_children: + max_bucket_size = 20 + if max_bucket_size >= 255: + raise ValueError("'max_bucket_size' must be smaller than 255") + + return max_children, max_bucket_size, bool(allow_duplicates) + + def _setup_c_archive(self): + # ffi.gc() ensures that the pointer is registered with the GC and the C + # destructor called. + return self._lib_archive_new( + self.dim, + self.max_children, + self.max_bucket_size, + self.allow_duplicates, + ) + + @staticmethod + def _c_get_x_values(c_archive) -> list[Any]: + if not lib.ndtree_has_x_values(c_archive): + return [] + x_list = ffi.new("void *[]", lib.ndtree_total_size(c_archive)) + lib.ndtree_get_x_values(c_archive, x_list) + x_list = [ffi.from_handle(x) for x in x_list] + return x_list + + def _lib_archive_insert_displaced(self, displaced) -> tuple[Any, bool]: + # lib.ndtree_insert_displaced() will free the current displaced, so + # Python does not need to do it. + ffi.gc(displaced, None) + self.displaced_p[0] = displaced + lib.ndtree_insert_displaced( + self._c_archive, self.displaced_p, self.all_displaced_p + ) + assert displaced != self.displaced_p[0] + displaced = self.displaced_p[0] + if displaced != ffi.NULL: + # lib.ndtree_insert_displaced() may create new memory, so we need ffi.gc(). + displaced = ffi.gc(displaced, self._lib_displaced_free) + + return displaced, self.all_displaced_p[0] + + +class TreapArchive(BaseUnboundedArchive): + # C functions + _lib_archive_add_and_get_displaced = lib.treap_archive_add + _lib_archive_dominated_by = lib.treap_archive_dominated_by + _lib_archive_dominates = lib.treap_archive_dominates + _lib_archive_find_exact_vector = lib.treap_archive_find_exact_vector + _lib_archive_find_x_values_of_vector = ( + lib.treap_archive_find_solutions_of_vector + ) + _lib_archive_free = lib.treap_archive_free + _lib_archive_get_contents = lib.treap_archive_get_contents + _lib_archive_get_unique_vectors = lib.treap_archive_get_unique_vectors + _lib_archive_has_x_values = lib.treap_archive_has_x_values + _lib_archive_new = lib.treap_archive_new + _lib_archive_total_size = lib.treap_archive_total_size + _lib_archive_unique_size = lib.treap_archive_unique_size + _lib_archive_iter_new = lib.treap_archive_iter_new + _lib_archive_iter_free = lib.treap_archive_iter_free + _lib_archive_iter_next = lib.treap_archive_iter_next + _lib_displaced_free = lib.treap_archive_free + _lib_displaced_get_unique_vectors = lib.treap_archive_get_unique_vectors + _lib_displaced_get_x_values = lib.treap_archive_get_x_values + _lib_displaced_total_len = lib.treap_archive_total_size + _lib_displaced_unique_len = lib.treap_archive_unique_size + + def __init__(self, maximise: bool | Sequence[bool]): + super().__init__(dim=2, maximise=maximise) + # We cache these pointers here so we don't keep allocating and freeing memory. + self.displaced_p = ffi.new("TreapArchive **") + + def _setup_c_archive(self): + return self._lib_archive_new() + + @staticmethod + def _c_get_x_values(c_archive) -> list[Any]: + if not lib.treap_archive_has_x_values(c_archive): + return [] + x_list = ffi.new("void *[]", lib.treap_archive_total_size(c_archive)) + lib.treap_archive_get_x_values(c_archive, x_list) + x_list = [ffi.from_handle(x) for x in x_list] + return x_list + + def _lib_archive_insert_displaced(self, displaced) -> tuple[Any, bool]: + # lib.treap_archive_insert_displaced() does not create new memory, so no ffi.gc(). + lib.treap_archive_insert_displaced( + self._c_archive, displaced, self.all_displaced_p + ) + return displaced, self.all_displaced_p[0] + + +# FIXME: implement maximise +@DocSubstitute() +def UnboundedArchive( # noqa: N802 + dim: int, + *, + maximise: bool | Sequence[bool] = False, + max_children: int = 0, + max_bucket_size: int = 0, + allow_duplicates: bool = True, +) -> BaseUnboundedArchive: + """Create an empty unbounded-size nondominated archive. The type of archive built depends on dim. + + A nondominated archive is a data structure that stores only nondominated points and allows efficient updates with new points. When a new point is added, if the new point is (weakly)-dominated by an existing point in the archive, the point would be rejected. If the new point is not (weakly)-dominated, it will be accepted. Points in the archive (weakly)-dominated by the new one will be removed. If the archive allows duplicates, then duplicated points, i.e., weakly-dominated but not dominated, are allowed in the archive. + + These archives are unbounded, that is, there is no limit in the number of nondominated points that can be stored in the archive. + + The archive can optionally store an arbitrary Python object associated to each point. When the archive allows duplicates, duplicated points may be associated to different Python objects. + + When :code:`dim = 2`, the archive is based on a treap data structure :footcite:p:`AraSei1989rst`. + When :code:`dim > 2`, the archive is based on a ND-Tree :footcite:p:`JasLus2018ndtree`. + + Parameters + ---------- + dim : + ${dim} + maximise : + ${maximise} + max_children : + Maximum number of children used in the ND-Tree (only when :code:`dim > 2`). + max_bucket_size : + Maximum number of solutions stored in each leaf of the ND-Tree (only when :code:`dim > 2`). + allow_duplicates : + Whether to allow or reject duplicated input points (i.e., objective vectors). + + Returns + ------- + Abstract base class that documents the methods available. + + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> arch = UnboundedArchive(dim=2) + >>> arch.add_many([[10, 10], [5, 20], [20, 5]]) + [True, True, True] + >>> arch.add([11, 11]) + False + + The ``len()`` function applied to the archive returns the number of points: + + >>> len(arch) + 3 + >>> arch.add([5, 5]) + True + >>> len(arch) + 1 + >>> arch.add_many([[6, 4], [4, 6]]) + [True, True] + + We can also iterate over the contents: + + >>> for z, _ in arch: + ... print(z) + [4. 6.] + [5. 5.] + [6. 4.] + + + It is possible to associate `x`-values to each point. + + >>> arch = UnboundedArchive(3) + >>> arch.add_many([[10, 10, 10], [5, 5, 20], [20, 5, 5]], x=["a", "b", "c"]) + [True, True, True] + >>> len(arch) + 3 + >>> arch.add([9, 9, 9], x="d") + True + >>> arch + [[20. 5. 5.] + [ 5. 5. 20.] + [ 9. 9. 9.]] + x = ['c', 'b', 'd'] + + The default is to allow duplicates: + + >>> arch.add([9, 9, 9], x="e") + True + >>> arch + [[20. 5. 5.] + [ 5. 5. 20.] + [ 9. 9. 9.] + [ 9. 9. 9.]] + x = ['c', 'b', 'd', 'e'] + + + The ``in`` keyword can be used to test if a point is in the archive: + + >>> [9, 9, 9] in arch + True + >>> [10, 10, 10] in arch + False + + And we can retrieve the object associated to a particular point: + + >>> arch.get_x_values_of_vector([9, 9, 9]) + ['d', 'e'] + >>> arch.get_x_values_of_vector([10, 10, 10]) + [] + + The contents of the archive can be obtained in different ways: + + >>> arch.get_contents() + (array([[20., 5., 5.], + [ 5., 5., 20.], + [ 9., 9., 9.], + [ 9., 9., 9.]]), ['c', 'b', 'd', 'e']) + >>> arch.get_unique_vectors() + array([[20., 5., 5.], + [ 5., 5., 20.], + [ 9., 9., 9.]]) + >>> arch.get_x_values() + ['c', 'b', 'd', 'e'] + >>> for z, x in arch: + ... print(f"{x} -> {z}") + c -> [20. 5. 5.] + b -> [ 5. 5. 20.] + d -> [9. 9. 9.] + e -> [9. 9. 9.] + + The vectors returned by iterating are read-only: + + >>> for z, x in arch: # doctest: +IGNORE_EXCEPTION_DETAIL +ELLIPSIS + ... z[0] = 1 + Traceback (most recent call last): + ValueError: assignment destination is read-only + + We can also mix maximisation and minimisation objectives: + + >>> arch = UnboundedArchive(dim=2, maximise=[False, True]) + >>> arch.add_many([[10, 10], [9, 9], [0, 0]]) + [True, True, True] + >>> arch.add([10, 9]) + False + >>> arch.dominated_by([0, 11]) + True + >>> arch.dominates([10, 11]) + False + >>> arch + [[ 0. 0.] + [ 9. 9.] + [10. 10.]] + x = [] + + """ + if dim > 2: + return NDTreeArchive( + dim, + maximise=maximise, + max_children=max_children, + max_bucket_size=max_bucket_size, + allow_duplicates=allow_duplicates, + ) + if dim == 2: + return TreapArchive(maximise=maximise) + if dim < 2: + raise ValueError("'dim' cannot be smaller than 2") + + raise NotImplementedError("This type of archive is not implemented yet") diff --git a/python/src/moocore/_docsubstitute.py b/python/src/moocore/_docsubstitute.py index 4e1f142d6..1a141b2ce 100644 --- a/python/src/moocore/_docsubstitute.py +++ b/python/src/moocore/_docsubstitute.py @@ -27,6 +27,8 @@ # ------------------------------------------------------------------------- "whv_rectangles": """Weighted rectangles that will bias the computation of the hypervolume. Maybe generated by :func:`eafdiff()` with ``rectangles=True`` or by :func:`choose_eafdiff()`.""", + # + "dim": """Dimension (number of criteria/objectives) of the input points stored in the archive.""", } diff --git a/python/src/moocore/_ffi_build.py b/python/src/moocore/_ffi_build.py index 87e32fff0..d907785c0 100644 --- a/python/src/moocore/_ffi_build.py +++ b/python/src/moocore/_ffi_build.py @@ -30,6 +30,8 @@ #include "whv.h" #include "whv_hype.h" #include "hvapprox.h" +#include "treap_archive.h" +#include "nd_tree.h" """ sources = [ "avl.c", @@ -49,6 +51,8 @@ "nondominated.c", "pareto.c", "rng.c", + "treap_archive.c", + "treap2d_fixed_size.c", "whv.c", "whv_hype.c", ] diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index ed3c932ad..6ffb0ae2e 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -2,9 +2,9 @@ import os from io import StringIO -from collections.abc import Callable, Sequence +from collections.abc import Callable from numpy.typing import ArrayLike, NDArray # For type hints -from typing import Literal, Any, NamedTuple, TypeAlias +from typing import Literal, Any, NamedTuple import numbers # NOTE: if we ever start using SciPy, we can use @@ -25,14 +25,14 @@ array_1d_of_length_n, is_integer_value, _get_seed_for_c, + _parse_maximise, + MaximiseType, ) from ._docsubstitute import DocSubstitute ## The CFFI library is used to create C bindings. from ._libmoocore import lib, ffi # ty: ignore[unresolved-import] -MaximiseType: TypeAlias = bool | Sequence[bool | int] | NDArray[np.bool_] - # Maximum number of objectives (columns) supported by the C library. These # mirror the limits defined in c/config.h: most functions support up to 255 # (the range of dimension_t == uint_fast8_t), while the hypervolume and @@ -171,11 +171,6 @@ def read_datasets(filename: str | os.PathLike[str] | StringIO) -> np.ndarray: return np.frombuffer(ffi.buffer(data_p, datasize)).reshape(-1, ncols_p[0]) -def _parse_maximise(maximise: MaximiseType, nobj: int) -> NDArray[np.bool_]: - """Convert maximise array or single bool to ndarray format.""" - return array_1d_of_length_n(maximise, nobj, name="maximise").astype(bool) - - def _parse_maximise_to_bool_array( maximise: MaximiseType, nobj: int ) -> ffi.CData: diff --git a/python/src/moocore/_utils.py b/python/src/moocore/_utils.py index 6ef8006e7..69c1cb999 100644 --- a/python/src/moocore/_utils.py +++ b/python/src/moocore/_utils.py @@ -1,9 +1,12 @@ -from numpy.typing import ArrayLike # For type hints -from typing import Any +from collections.abc import Sequence +from numpy.typing import ArrayLike, NDArray # For type hints +from typing import Any, TypeAlias import numpy as np from ._libmoocore import ffi # ty: ignore[unresolved-import] +MaximiseType: TypeAlias = bool | Sequence[bool | int] | NDArray[np.bool_] + def get1_and_is_copied(x: Any, x_: Any) -> tuple[Any, bool]: x_copied = id(x) != id(x_) @@ -27,6 +30,32 @@ def unique_nosort(array: ArrayLike, axis: int | None = None) -> np.ndarray: return uniq[index.argsort()] +def _as_1d_float_array( + f: ArrayLike, dim: int | None = None, name: str = "f" +) -> np.ndarray: + """Convert input to a validated 1D float array.""" + result = np.ascontiguousarray(f, dtype=float) + if result.ndim != 1: + raise ValueError(f"'{name}' must be a 1D objective vector") + if dim is not None and result.shape[0] != dim: + raise ValueError(f"'{name}' must have length {dim}") + return result + + +def _as_2d_float_array( + f: ArrayLike, *, dim: int | None = None, name: str = "f" +) -> np.ndarray: + """Convert input to a validated 2D float matrix.""" + result = np.ascontiguousarray(f, dtype=float) + if result.ndim == 1: + result = result.reshape(1, -1) + if result.ndim != 2: + raise ValueError(f"'{name}' must be a 1D vector or 2D matrix") + if dim is not None and result.shape[1] != dim: + raise ValueError(f"'{name}' must have {dim} columns") + return result + + def np2d_to_double_array( x: np.ndarray, ctype_shape: tuple[str, str] = ("int", "int") ) -> tuple[ffi.CData, ffi.CData, ffi.CData]: @@ -87,3 +116,8 @@ def _get_seed_for_c(seed: int | np.random.Generator | None) -> ffi.CData: if not is_integer_value(seed): seed = np.random.default_rng(seed).integers(2**32 - 2, dtype=np.uint32) return ffi.cast("uint32_t", seed) + + +def _parse_maximise(maximise: MaximiseType, nobj: int) -> NDArray[np.bool_]: + """Convert maximise array or single bool to ndarray format.""" + return array_1d_of_length_n(maximise, nobj, name="maximise").astype(bool) diff --git a/python/src/moocore/libmoocore.h b/python/src/moocore/libmoocore.h index 5f0c0a45e..253b839ff 100644 --- a/python/src/moocore/libmoocore.h +++ b/python/src/moocore/libmoocore.h @@ -77,6 +77,83 @@ double hv_approx_rphi_fang_wang_plus(const double * restrict data, const double * restrict ref, const boolvec * restrict maximise, uint_fast32_t nsamples); + +// archive.h +typedef enum archive_insert_result_t { + ARCHIVE_INSERT_X_MUST_BE_NOT_NULL =..., + ARCHIVE_INSERT_X_MUST_BE_NULL, + ARCHIVE_INSERT_MEMORY_ERROR, + ARCHIVE_INSERT_REJECTED, + ARCHIVE_INSERT_ACCEPTED, + ARCHIVE_INSERT_DUPLICATED, +}; +typedef struct VoidList { + uint32_t n; + uint32_t cap; + const void ** list; +} VoidList; + +// treap_archive.h +typedef struct TreapArchive TreapArchive; +typedef struct TreapNode TreapNode; + +TreapArchive * treap_archive_new(void); +void treap_archive_free(TreapArchive * t); +bool treap_archive_find_exact_vector(const TreapArchive * t, const double * z); +VoidList * treap_archive_find_solutions_of_vector(const TreapArchive * restrict tree, const double * restrict z); +bool treap_archive_dominated_by(const TreapArchive * self, const double * z); +bool treap_archive_dominates(const TreapArchive * self, const double * z); +int treap_archive_add(TreapArchive *self, const double * z, const void * x, bool check, TreapArchive **displaced); +void treap_archive_insert_displaced(TreapArchive *t, TreapArchive *incoming, bool *all_displaced); +void treap_archive_get_unique_vectors(const TreapArchive *t, const double ** z_list); +void treap_archive_get_x_values(const TreapArchive *t, const void ** x_list); +bool treap_archive_has_x_values(const TreapArchive *t); +void treap_archive_get_contents(const TreapArchive *t, double * z_list, const void ** x_list); +size_t treap_archive_total_size(const TreapArchive *t); +size_t treap_archive_unique_size(const TreapArchive *t); + +typedef struct TreapIterator TreapIterator; +TreapIterator *treap_archive_iter_new(const TreapArchive *); +int treap_archive_iter_next(TreapIterator *it, const double **z, const void **x); +void treap_archive_iter_free(TreapIterator *it); + +// nd_tree.h +typedef struct FlexBucket FlexBucket; +void flex_bucket_free(FlexBucket * bucket); +size_t flex_bucket_unique_len(const FlexBucket * bucket); +size_t flex_bucket_total_len(const FlexBucket * bucket); +void flex_bucket_get_unique_vectors(const FlexBucket * bucket, const double ** z_list); +void flex_bucket_get_x_values(const FlexBucket * bucket, const void ** x_list); + +typedef struct NDTreeConfig { + uint8_t max_children; + uint8_t max_bucket_size; + bool allow_duplicates; +} NDTreeConfig; + +typedef struct NDTreeArchive NDTreeArchive; +typedef struct NDTreeNode NDTreeNode; +NDTreeArchive * ndtree_new(dimension_t dim, uint8_t max_children, uint8_t max_bucket_size, bool allow_duplicates); +void ndtree_free(NDTreeArchive *self); +bool ndtree_has_x_values(const NDTreeArchive *t); +size_t ndtree_total_size(const NDTreeArchive *self); +size_t ndtree_unique_size(const NDTreeArchive *self); +bool ndtree_dominates(const NDTreeArchive *self, const double *z); +bool ndtree_dominated_by(const NDTreeArchive * self, const double * z); +int ndtree_add(NDTreeArchive *self, const double *z, const void *x, bool check, FlexBucket **displaced); +int ndtree_insert_displaced(NDTreeArchive *tree, FlexBucket **displaced_p, bool *all_displaced); +bool ndtree_find_exact_vector(const NDTreeArchive *self, const double *z); +VoidList * ndtree_find_solutions_of_vector(const NDTreeArchive * restrict tree, const double * restrict z); +void ndtree_get_contents(const NDTreeArchive *self, double * z_list, const void ** x_list); +void ndtree_get_unique_vectors(const NDTreeArchive *self, const double ** z_list); +void ndtree_get_x_values(const NDTreeArchive *arch, const void ** x_list); + +typedef struct NDTreeIterator NDTreeIterator; +NDTreeIterator *ndtree_iter_new(const NDTreeArchive *); +int ndtree_iter_next(NDTreeIterator *it, const double **z, const void **x); +void ndtree_iter_free(NDTreeIterator *it); + + /* typedef ... hype_sample_dist; hype_sample_dist * hype_dist_unif_new(unsigned long seed); diff --git a/python/tests/conftest.py b/python/tests/conftest.py index fd1ff8e1a..715513067 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -70,3 +70,21 @@ def _call(fn, *args, immutable=None, **kwargs): return result return _call + + +class AutoIncrement: + """Useful for counting.""" + + def __init__(self, start=0): + self._x = start - 1 + + def __call__(self): + """Return current value and increment counter.""" + self._x += 1 + return self._x + + +@pytest.fixture +def count(): + """Instantiate a counter of type AutoIncrement().""" + return AutoIncrement() diff --git a/python/tests/test_archiving.py b/python/tests/test_archiving.py new file mode 100644 index 000000000..2e554c2e7 --- /dev/null +++ b/python/tests/test_archiving.py @@ -0,0 +1,468 @@ +# ruff: noqa: D100, D103, N802 +import pytest +import numpy as np +from numpy.testing import ( + assert_array_equal, +) +import moocore + + +def no_print(x): + return + + +def extend_point(z, dim): + z = np.asarray(z) + z_ext = np.zeros(dim) + z_ext[-len(z) :] = z + return z_ext + + +def test_UnboundedArchive2D_add(): + archive = moocore.UnboundedArchive(dim=2) + assert archive.add([1, 10]) + assert len(archive) == 1 + assert archive.add([10, 1]) + assert len(archive) == 2 + assert archive.add([5, 5]) + assert len(archive) == 3 + assert archive.unique_len() == 3 + # Duplicates increase len() + assert archive.add([5, 5]) + assert len(archive) == 4 + assert archive.unique_len() == 3 + assert archive.add([4, 4]) + assert len(archive) == 3 + assert not archive.add([5, 5]) + assert len(archive) == 3 + assert archive.add([1, 1]) + assert len(archive) == 1 + + # Test mixing up x = None and x not None. + with pytest.raises(ValueError, match=r"'x' must be None"): + archive.add([1, 1], x=2) + + with pytest.raises(ValueError, match=r"'x' must be None"): + archive.add([0, 0], x=2) + + archive = moocore.UnboundedArchive(dim=2) + archive.add([1, 10], x="a") + with pytest.raises(ValueError, match=r"'x' cannot be None"): + archive.add([10, 1]) + + +def test_UnboundedArchive2D_dominance(): + archive = moocore.UnboundedArchive(dim=2) + assert_array_equal( + archive.add_many([[1, 10], [5, 5], [10, 1]]), [True, True, True] + ) + # Test __contains__ + assert [5, 5] in archive + assert [4, 4] not in archive + + assert not archive.dominated_by([4, 4]) + assert not archive.dominated_by([5, 5]) + assert not archive.dominated_by([6, 6]) + assert archive.dominated_by([1, 1]) + + assert archive.dominates([6, 6]) # One point in the archive dominates. + assert archive.dominates([10, 10]) # Whole archive dominates. + assert not archive.dominates( + [4, 4] + ) # One point in the archive is dominated. + assert not archive.dominates( + [5, 5] + ) # One point in the archive is weakly dominated. + assert not archive.dominates([6, 4]) # Mutually non-dominated + assert not archive.dominates([1, 1]) + + assert archive.add([1, 1]) + assert len(archive) == 1 + assert not archive.dominated_by( + [1, 1] + ) # Duplicated do not count as dominance. + assert not archive.dominated_by([2, 1]) + assert not archive.dominated_by([1, 2]) + assert archive.dominated_by([0, 1]) + assert archive.dominated_by([1, 0]) + + assert not archive.dominates( + [1, 1] + ) # Duplicated do not count as dominance. + assert not archive.dominates([1, 0]) + assert not archive.dominates([0, 1]) + assert archive.dominates([2, 1]) + assert archive.dominates([1, 2]) + + +def test_UnboundedArchive2D_displaced(count): + dim = 2 + pf = moocore.UnboundedArchive(dim) + + accepted, displaced = pf._add_with_displaced_raw( + [1, 10], x=count(), check=False, discard=True + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [10, 1], x=count(), check=False, discard=True + ) + assert accepted and displaced is None + + assert pf.dominated_by([1, 1]) + assert not pf.dominated_by([2, 9]) + + accepted, displaced = pf._add_with_displaced_raw( + [2, 9], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + + assert not pf.dominated_by([2, 9]) + assert not pf.dominated_by([9, 2]) + assert not pf.dominates([2, 9]) + assert not pf.dominates([9, 2]) + assert pf.dominates([3, 10]) + + accepted, displaced = pf._add_with_displaced_raw( + [3, 8], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [3, 8], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [2, 10], x=count(), check=True, discard=False + ) + assert not accepted and displaced is None + accepted, displaced_1 = pf._add_with_displaced_raw( + [0, 10], x=count(), check=True, discard=False + ) + assert accepted and displaced_1 is not None + accepted, displaced_2 = pf._add_with_displaced_raw( + [2, 8], x=count(), check=True, discard=False + ) + assert accepted and displaced_2 is not None + + z, x = pf.get_contents() + no_print(z) + no_print(x) + + pf2 = moocore.UnboundedArchive(dim) + accepted, displaced = pf2._add_with_displaced_raw( + [0, 11], x=count(), check=False, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [11, 1], x=count(), check=False, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [10, 3], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [4, 9], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + z, x = pf2.get_contents() + no_print(z) + no_print(x) + displaced = pf2._add_displaced_raw(displaced_1) + assert displaced is None + displaced = pf2._add_displaced_raw(displaced_2) + assert displaced is not None + no_print(displaced) + z, x = pf2.get_contents() + no_print(z) + no_print(x) + + +@pytest.mark.parametrize("dim", [3, 4, 5, 10]) +def test_UnboundedArchive_add(dim): + + archive = moocore.UnboundedArchive(dim=dim) + + def archive_add_extended(z, **kwargs): + z_dim = extend_point(z, dim) + return archive.add(z_dim, **kwargs) + + def assert_add(z, expect_res, expect_len): + assert archive_add_extended(z) == expect_res + assert len(archive) == expect_len + + assert_add([1, 10, 1], True, 1) + assert_add([10, 1, 10], True, 2) + assert_add([5, 5, 5], True, 3) + assert archive.unique_len() == 3 + # Duplicates increase len() + assert_add([5, 5, 5], True, 4) + assert archive.unique_len() == 3 + assert_add([4, 4, 4], True, 3) + assert_add([5, 5, 5], False, 3) + assert_add([1, 1, 1], True, 1) + + # Test mixing up x = None and x not None. + with pytest.raises(ValueError, match=r"'x' must be None"): + archive_add_extended([2, 0, 1], x=2) + + with pytest.raises(ValueError, match=r"'x' must be None"): + archive_add_extended([1, 1, 1], x=2) + + with pytest.raises(ValueError, match=r"'x' must be None"): + archive_add_extended([0, 0, 0], x=2) + + archive = moocore.UnboundedArchive(dim=dim) + archive_add_extended([1, 10, 1], x="a") + with pytest.raises(ValueError, match=r"'x' cannot be None"): + archive_add_extended([10, 1, 10]) + + +# FIXME: How to generalize this test to any number of dimensions? +def test_UnboundedArchive3D_dominance(): + + archive = moocore.UnboundedArchive(dim=3) + assert_array_equal( + archive.add_many([[1, 10, 1], [5, 5, 5], [10, 1, 10]]), + [True, True, True], + ) + # Test __contains__ + assert [5, 5, 5] in archive + assert [4, 4, 4] not in archive + + assert not archive.dominated_by([4, 4, 4]) + assert not archive.dominated_by([5, 5, 5]) + assert not archive.dominated_by([6, 6, 6]) + assert archive.dominated_by([1, 1, 1]) + + assert archive.dominates([6, 6, 6]) # One point in the archive dominates. + assert archive.dominates([10, 10, 10]) # Whole archive dominates. + assert not archive.dominates( + [4, 4, 4] + ) # One point in the archive is dominated. + assert not archive.dominates( + [5, 5, 5] + ) # One point in the archive is weakly dominated. + assert not archive.dominates([11, 1, 1]) # Mutually non-dominated + assert not archive.dominates([1, 1, 1]) + + assert archive.add([1, 1, 1]) + assert len(archive) == 1 + assert not archive.dominated_by( + [1, 1, 1] + ) # Duplicated do not count as dominance. + assert not archive.dominated_by([2, 1, 1]) + assert not archive.dominated_by([1, 2, 1]) + assert archive.dominated_by([0, 1, 0]) + assert archive.dominated_by([1, 0, 0]) + + assert not archive.dominates( + [1, 1, 1] + ) # Duplicated do not count as dominance. + assert not archive.dominates([1, 0, 0]) + assert not archive.dominates([0, 1, 0]) + assert archive.dominates([2, 1, 20]) + assert archive.dominates([1, 2, 4]) + + +# FIXME: How to generalize this test to any number of dimensions? +def test_UnboundedArchive3D_displaced(count): + dim = 3 + pf = moocore.UnboundedArchive(dim) + + accepted, displaced = pf._add_with_displaced_raw( + [1, 11, 10], x=count(), check=False, discard=True + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [10, 1, 10], x=count(), check=False, discard=True + ) + assert accepted and displaced is None + + assert pf.dominated_by([1, 1, 1]) + assert not pf.dominated_by([2, 9, 1]) + + accepted, displaced = pf._add_with_displaced_raw( + [2, 9, 1], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + + assert not pf.dominated_by([2, 9, 1]) + assert not pf.dominated_by([9, 2, 1]) + assert not pf.dominates([2, 9, 1]) + assert not pf.dominates([9, 2, 1]) + assert pf.dominates([3, 10, 10]) + + accepted, displaced = pf._add_with_displaced_raw( + [3, 8, 10], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [3, 8, 10], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf._add_with_displaced_raw( + [2, 10, 10], x=count(), check=True, discard=False + ) + assert not accepted and displaced is None + accepted, displaced_1 = pf._add_with_displaced_raw( + [0, 10, 9], x=count(), check=True, discard=False + ) + assert accepted and displaced_1 is not None + accepted, displaced_2 = pf._add_with_displaced_raw( + [2, 8, 1], x=count(), check=True, discard=False + ) + assert accepted and displaced_2 is not None + + z, x = pf.get_contents() + no_print(z) + no_print(x) + + pf2 = moocore.UnboundedArchive(dim) + accepted, displaced = pf2._add_with_displaced_raw( + [ + 1, + 10, + 11, + ], + x=count(), + check=False, + discard=False, + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [10, 1, 11], x=count(), check=False, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [11, 3, 9], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + accepted, displaced = pf2._add_with_displaced_raw( + [4, 11, 8], x=count(), check=True, discard=False + ) + assert accepted and displaced is None + z, x = pf2.get_contents() + no_print(z) + no_print(x) + displaced = pf2._add_displaced_raw(displaced_1) + assert displaced is None + displaced = pf2._add_displaced_raw(displaced_2) + assert displaced is not None + no_print(displaced) + z, x = pf2.get_contents() + no_print(z) + no_print(x) + + +@pytest.mark.parametrize("dim", [2, 3, 5, 10]) +def test_UnboundedArchive_random(dim, count): + """Randomized tests for any number of dimensions.""" + seed = np.random.default_rng().integers(2**32 - 2) + print(f"seed={seed}") + rng = np.random.default_rng(seed) + points = rng.integers(0, 100, (1000, dim)) + archive = moocore.UnboundedArchive(dim) + for p in points: + archive.add(p, x=count()) + + z, x = archive.get_contents() + points = moocore.filter_dominated(points, keep_weakly=True) + z = z.astype(int) + assert_array_equal(points, z[np.argsort(x), :]) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_invalid_iterator(dim): + archive = moocore.UnboundedArchive(dim) + rng = np.random.default_rng() + z = rng.uniform(size=(3, dim)) + archive.add_many(z) + + with pytest.raises( + RuntimeError, match=r"invalid 'BaseUnboundedArchive' iterator" + ): + for z, x in archive: + archive.add(np.full(dim, 0)) + + # archive.add() succeeded. + assert len(archive) == 1 + + # Test that the iterator keeps the archive alive. + iterator = iter(archive) + del archive + z, x = next(iterator) + assert len(z) == dim + assert x is None + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_maximise(dim): + maximise = np.array([True, False, True]) + n = 10 + z = np.random.uniform(size=(n, dim)) + z_orig = z.copy() + + arch = moocore.UnboundedArchive(dim=dim, maximise=maximise[:dim]) + arch.add_many(z, x=np.arange(len(z))) + for z, x in arch: + assert np.array_equal(z, z_orig[x, :]) + for z, x in arch: + assert np.array_equal(z, z_orig[x, :]) + + z, x = arch.get_contents() + assert np.array_equal(z, z_orig[np.array(x), :]) + + z, x = arch.get_contents() + assert np.array_equal(z, z_orig[np.array(x), :]) + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_find_exact(dim): + z = moocore.generate_ndset(1000, dim, method="simplex", integer=True) + arch = moocore.UnboundedArchive(dim=dim) + arch.add_many(z, x=np.arange(len(z))) + for i in [1, 500, 990]: + res = arch.get_x_values_of_vector(z[i, :]) + assert res[0] == i + + y = z.max(axis=0) + res = arch.get_x_values_of_vector(y) + assert len(res) == 0 + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_empty_archive(dim): + """Test that operations on an empty archive works as expected.""" + arch = moocore.UnboundedArchive(dim=dim) + z = np.full(dim, 1) + + assert len(arch) == 0 + assert arch.unique_len() == 0 + assert not arch.has_x_values() + assert z not in arch + assert not arch.dominates(z) + assert not arch.dominated_by(z) + x = arch.get_x_values_of_vector(z) + assert len(x) == 0 + for z, _ in arch: + assert False + z, x = arch.get_contents() + assert len(z) == 0 + assert len(x) == 0 + z = arch.get_unique_vectors() + assert len(z) == 0 + x = arch.get_x_values() + assert len(x) == 0 + + +def test_ndtree_dominated_by(): + dim = 3 + arch = moocore.UnboundedArchive(dim=dim) + n = 1000 + z = np.random.uniform(size=(n, dim)) + # We cannot use add_many() because it may filter before insertion. + for z_i in z: + arch.add(z_i) + nondom = moocore.filter_dominated(z) + ideal = nondom.min(axis=0) + assert arch.dominated_by(ideal) diff --git a/r/inst/REFERENCES.bib b/r/inst/REFERENCES.bib index b752aab29..871755195 100644 --- a/r/inst/REFERENCES.bib +++ b/r/inst/REFERENCES.bib @@ -298,6 +298,42 @@ @article{HerWer1987tabucol pages = {345--351} } +@article{JasLus2018ndtree, + author = { Andrzej Jaszkiewicz and Thibaut Lust }, + title = {{ND}-Tree-Based Update: A Fast Algorithm for the Dynamic + Nondominance Problem}, + journal = {IEEE Transactions on Evolutionary Computation}, + year = 2018, + volume = 22, + number = 5, + pages = {778--791}, + doi = {10.1109/TEVC.2018.2799684}, + publisher = {IEEE}, + abstract = {In this paper, we propose a new method called ND-Tree-based + update (ND-Tree) for the dynamic nondominance problem, i.e., + the problem of online update of a Pareto archive composed of + mutually nondominated points. It uses a new ND-Tree data + structure in which each node represents a subset of points + contained in a hyperrectangle defined by its local + approximate ideal and nadir points. By building subsets + containing points located close in the objective space and + using basic properties of the local ideal and nadir points we + can efficiently avoid searching many branches in the + tree. ND-Tree may be used in multiobjective evolutionary + algorithms and other multiobjective metaheuristics to update + an archive of potentially nondominated points. We prove that + the proposed algorithm has sublinear time complexity under + mild assumptions. We experimentally compare ND-Tree to the + simple list, Quad-tree, and M-Front methods using artificial + and realistic benchmarks with up to ten objectives and show + that with this new method substantial reduction of the number + of point comparisons and computational time can be + obtained. Furthermore, we apply the method to the + nondominated sorting problem showing that it is highly + competitive to some recently proposed algorithms dedicated to + this problem.} +} + @article{Jen03, author = {M. T. Jensen}, title = {Reducing the run-time complexity of multiobjective {EA}s: The @@ -455,6 +491,16 @@ @article{ZitThiLauFon2003:tec doi = {10.1109/TEVC.2003.810758} } +@inproceedings{AraSei1989rst, + author = {Aragon, C. R. and Seidel, R. G.}, + title = {Randomized search trees}, + booktitle = {30th Annual Symposium on Foundations of Computer Science}, + year = 1989, + pages = {540--545}, + keywords = {Binary search trees, Treap, Dynamic programming}, + doi = {10.1109/SFCS.1989.63531} +} + @incollection{AugBadBroZit2009gecco, booktitle = {Proceedings of the Genetic and Evolutionary Computation Conference, GECCO 2009}, address = { New York, NY}, diff --git a/r/src/Makevars b/r/src/Makevars index 8164b3bdd..9a6ee258b 100644 --- a/r/src/Makevars +++ b/r/src/Makevars @@ -24,6 +24,7 @@ MOOCORE_SRC_FILES = \ pareto.c \ r2_exact.c \ rng.c \ + treap2d_fixed_size.c \ whv.c \ whv_hype.c diff --git a/r/src/install.libs.R b/r/src/install.libs.R index ac4a6e39f..950ba3717 100644 --- a/r/src/install.libs.R +++ b/r/src/install.libs.R @@ -1,4 +1,5 @@ files <- c( + "archiver", "dominatedsets", "eaf", "epsilon",