From ee9a8b1f0b13ba5a1ccc030168ff092f301e6922 Mon Sep 17 00:00:00 2001 From: Jakob Bossek Date: Wed, 13 May 2015 15:52:11 +0200 Subject: [PATCH] add function to compute dominated hypervolume (not finished) --- DESCRIPTION | 1 + NAMESPACE | 2 + R/computeDominatedHypervolume.R | 37 ++ R/zzz.R | 3 +- man/computeDominatedHypervolume.Rd | 25 ++ src/Makevars | 7 + src/avl.c | 589 +++++++++++++++++++++++++++++ src/avl.h | 186 +++++++++ src/avl.o | Bin 0 -> 21276 bytes src/computeHypervolume.c | 29 ++ src/computeHypervolume.o | Bin 0 -> 3440 bytes src/ecr.so | Bin 0 -> 25428 bytes src/hv.c | 525 +++++++++++++++++++++++++ src/hv.h | 48 +++ src/hv.o | Bin 0 -> 17320 bytes src/macros.h | 22 ++ 16 files changed, 1473 insertions(+), 1 deletion(-) create mode 100644 R/computeDominatedHypervolume.R create mode 100644 man/computeDominatedHypervolume.Rd create mode 100644 src/Makevars create mode 100644 src/avl.c create mode 100644 src/avl.h create mode 100644 src/avl.o create mode 100644 src/computeHypervolume.c create mode 100644 src/computeHypervolume.o create mode 100755 src/ecr.so create mode 100644 src/hv.c create mode 100644 src/hv.h create mode 100644 src/hv.o create mode 100644 src/macros.h diff --git a/DESCRIPTION b/DESCRIPTION index ff78134..7842b3f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,3 +27,4 @@ Imports: Suggests: testthat (>= 0.9.1) ByteCompile: yes +LazyData: yes diff --git a/NAMESPACE b/NAMESPACE index 8f1632d..3b740e3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,7 @@ S3method(print,ecr_control) S3method(print,ecr_result) export(computeAverageHausdorffDistance) export(computeCrowdingDistance) +export(computeDominatedHypervolume) export(computeGenerationalDistance) export(computeInvertedGenerationalDistance) export(doNondominatedSorting) @@ -63,3 +64,4 @@ import(gridExtra) import(parallelMap) import(reshape2) import(smoof) +useDynLib(ecr) diff --git a/R/computeDominatedHypervolume.R b/R/computeDominatedHypervolume.R new file mode 100644 index 0000000..340148c --- /dev/null +++ b/R/computeDominatedHypervolume.R @@ -0,0 +1,37 @@ +#' Computation of the dominated hypervolume. +#' +#' Given a set of points \code{x}, this function computes the dominated hypervolume +#' of the points regarding the reference point \code{ref.point}. If the latter +#' is not provided, one is automatically determined by computing the maximum +#' in each dimension. +#' +#' @param x [\code{matrix}]\cr +#' Matrix of points. +#' @param ref.point [\code{numeric} | \code{NULL}]\cr +#' Reference point. Set to the maximum in each dimension by default if not provided. +#' @return [\code{numeric(1)}] Dominated hypervolume. +#' @export +computeDominatedHypervolume = function(x, ref.point = NULL) { + # sanity checks + assertMatrix(x, mode = "numeric", any.missing = FALSE) + if (any(is.infinite(x))) { + warningf("Set of points contains infinite %i values.", which(is.infinite(x))) + return(NaN) + } + + if (is.null(ref.point)) { + ref.point = apply(x, 1, max) + } + + if (length(ref.point) != nrow(x)) { + stopf("Set of points and reference point need to have the same dimension, but + set of points has dimension %i and reference points has dimension %i.", nrow(x), length(ref.point)) + } + + if (any(is.infinite(ref.point))) { + warningf("Reference point contains infinite %i values.", which(is.infinite(ref.point))) + return(NaN) + } + + return(.Call("computeDominatedHypervolumeC", as.numeric(x), ref.point)) +} diff --git a/R/zzz.R b/R/zzz.R index 2fb07fb..3c57d2a 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -6,4 +6,5 @@ #' @import parallelMap #' @import reshape2 #' @import gridExtra -NULL \ No newline at end of file +#' @useDynLib ecr +NULL diff --git a/man/computeDominatedHypervolume.Rd b/man/computeDominatedHypervolume.Rd new file mode 100644 index 0000000..4f8d0ba --- /dev/null +++ b/man/computeDominatedHypervolume.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/computeDominatedHypervolume.R +\name{computeDominatedHypervolume} +\alias{computeDominatedHypervolume} +\title{Computation of the dominated hypervolume.} +\usage{ +computeDominatedHypervolume(x, ref.point = NULL) +} +\arguments{ +\item{x}{[\code{matrix}]\cr +Matrix of points.} + +\item{ref.point}{[\code{numeric} | \code{NULL}]\cr +Reference point. Set to the maximum in each dimension by default if not provided.} +} +\value{ +[\code{numeric(1)}] Dominated hypervolume. +} +\description{ +Given a set of points \code{x}, this function computes the dominated hypervolume +of the points regarding the reference point \code{ref.point}. If the latter +is not provided, one is automatically determined by computing the maximum +in each dimension. +} + diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..66950e6 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,7 @@ +all: $(SHLIB) + +hv.o: hv.c + $(CC) $(ALL_CPPFLAGS) $(ALL_CFLAGS) -DVARIANT=4 -c -o hv.o hv.c + +avl.o: avl.c + $(CC) $(ALL_CPPFLAGS) $(ALL_CFLAGS) -DVARIANT=4 -c -o avl.o avl.c diff --git a/src/avl.c b/src/avl.c new file mode 100644 index 0000000..5771b97 --- /dev/null +++ b/src/avl.c @@ -0,0 +1,589 @@ +/***************************************************************************** + + avl.c - 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. + +*****************************************************************************/ + +#include +#include +#include +#include "avl.h" + +static void avl_rebalance(avl_tree_t *, avl_node_t *); + +#ifdef AVL_COUNT +#define NODE_COUNT(n) ((n) ? (n)->count : 0) +#define L_COUNT(n) (NODE_COUNT((n)->left)) +#define R_COUNT(n) (NODE_COUNT((n)->right)) +#define CALC_COUNT(n) (L_COUNT(n) + R_COUNT(n) + 1) +#endif + +#ifdef AVL_DEPTH +#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) ((L_DEPTH(n)>R_DEPTH(n)?L_DEPTH(n):R_DEPTH(n)) + 1) +#endif + +#ifndef AVL_DEPTH +/* Also known as ffs() (from BSD) */ +static int lg(unsigned int u) { + int r = 1; + if(!u) return 0; + if(u & 0xffff0000) { u >>= 16; r += 16; } + if(u & 0x0000ff00) { u >>= 8; r += 8; } + if(u & 0x000000f0) { u >>= 4; r += 4; } + if(u & 0x0000000c) { u >>= 2; r += 2; } + if(u & 0x00000002) r++; + return r; +} +#endif + +static int avl_check_balance(avl_node_t *avlnode) { +#ifdef AVL_DEPTH + int d; + d = R_DEPTH(avlnode) - L_DEPTH(avlnode); + return d<-1?-1:d>1?1:0; +#else +/* int d; + * d = lg(R_COUNT(avlnode)) - lg(L_COUNT(avlnode)); + * d = d<-1?-1:d>1?1:0; + */ +#ifdef AVL_COUNT + int pl, r; + + pl = lg(L_COUNT(avlnode)); + r = R_COUNT(avlnode); + + if(r>>pl+1) + return 1; + if(pl<2 || r>>pl-2) + return 0; + return -1; +#else +#error No balancing possible. +#endif +#endif +} + +#ifdef AVL_COUNT +unsigned int avl_count(const avl_tree_t *avltree) { + return NODE_COUNT(avltree->top); +} + +avl_node_t *avl_at(const avl_tree_t *avltree, unsigned int index) { + avl_node_t *avlnode; + unsigned int c; + + avlnode = avltree->top; + + while(avlnode) { + c = L_COUNT(avlnode); + + if(index < c) { + avlnode = avlnode->left; + } else if(index > c) { + avlnode = avlnode->right; + index -= c+1; + } else { + return avlnode; + } + } + return NULL; +} + +unsigned int avl_index(const avl_node_t *avlnode) { + avl_node_t *next; + unsigned int c; + + c = L_COUNT(avlnode); + + while((next = avlnode->parent)) { + if(avlnode == next->right) + c += L_COUNT(next) + 1; + avlnode = next; + } + + return c; +} +#endif + +int avl_search_closest(const avl_tree_t *avltree, const void *item, avl_node_t **avlnode) { + avl_node_t *node; + avl_compare_t cmp; + int c; + + if(!avlnode) + avlnode = &node; + + node = avltree->top; + + if(!node) + return *avlnode = NULL, 0; + + cmp = avltree->cmp; + + for(;;) { + 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; + } + } +} + +/* + * avl_search: + * Return a pointer to a node with the given item in the tree. + * If no such item is in the tree, then NULL is returned. + */ +avl_node_t *avl_search(const avl_tree_t *avltree, const void *item) { + avl_node_t *node; + return avl_search_closest(avltree, item, &node) ? NULL : node; +} + +avl_tree_t *avl_init_tree(avl_tree_t *rc, avl_compare_t cmp, avl_freeitem_t freeitem) { + if(rc) { + rc->head = NULL; + rc->tail = NULL; + rc->top = NULL; + rc->cmp = cmp; + rc->freeitem = freeitem; + } + return rc; +} + +avl_tree_t *avl_alloc_tree(avl_compare_t cmp, avl_freeitem_t freeitem) { + return avl_init_tree(malloc(sizeof(avl_tree_t)), cmp, freeitem); +} + +void avl_clear_tree(avl_tree_t *avltree) { + avltree->top = avltree->head = avltree->tail = NULL; +} + +void avl_free_nodes(avl_tree_t *avltree) { + avl_node_t *node, *next; + avl_freeitem_t freeitem; + + freeitem = avltree->freeitem; + + for(node = avltree->head; node; node = next) { + next = node->next; + if(freeitem) + freeitem(node->item); + free(node); + } + avl_clear_tree(avltree); +} + +/* + * avl_free_tree: + * Free all memory used by this tree. If freeitem is not NULL, then + * it is assumed to be a destructor for the items referenced in the avl_ + * tree, and they are deleted as well. + */ +void avl_free_tree(avl_tree_t *avltree) { + avl_free_nodes(avltree); + free(avltree); +} + +static void avl_clear_node(avl_node_t *newnode) { + newnode->left = newnode->right = NULL; + #ifdef AVL_COUNT + newnode->count = 1; + #endif + #ifdef AVL_DEPTH + newnode->depth = 1; + #endif +} + +avl_node_t *avl_init_node(avl_node_t *newnode, void *item) { + if(newnode) { + avl_clear_node(newnode); + newnode->item = item; + } + return newnode; +} + +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; +} + +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); + + 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; +} + +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); + + 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; +} + +avl_node_t *avl_insert_node(avl_tree_t *avltree, avl_node_t *newnode) { + avl_node_t *node; + + if(!avltree->top) + return avl_insert_top(avltree, newnode); + + switch(avl_search_closest(avltree, newnode->item, &node)) { + case -1: + return avl_insert_before(avltree, node, newnode); + case 1: + return avl_insert_after(avltree, node, newnode); + } + + return NULL; +} + +/* + * avl_insert: + * Create a new node and insert an item there. + * Returns the new node on success or NULL if no memory could be allocated. + */ +avl_node_t *avl_insert(avl_tree_t *avltree, void *item) { + avl_node_t *newnode; + + newnode = avl_init_node(malloc(sizeof(avl_node_t)), item); + if(newnode) { + if(avl_insert_node(avltree, newnode)) + return newnode; + free(newnode); + errno = EEXIST; + } + return NULL; +} + +/* + * avl_unlink_node: + * Removes the given node. 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.) + */ +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); +} + +void *avl_delete_node(avl_tree_t *avltree, avl_node_t *avlnode) { + void *item = NULL; + if(avlnode) { + item = avlnode->item; + avl_unlink_node(avltree, avlnode); + if(avltree->freeitem) + avltree->freeitem(item); + free(avlnode); + } + return item; +} + +void *avl_delete(avl_tree_t *avltree, const void *item) { + return avl_delete_node(avltree, avl_search(avltree, item)); +} + +avl_node_t *avl_fixup_node(avl_tree_t *avltree, avl_node_t *newnode) { + avl_node_t *oldnode = NULL, *node; + + if(!avltree || !newnode) + return NULL; + + node = newnode->prev; + if(node) { + oldnode = node->next; + node->next = newnode; + } else { + avltree->head = newnode; + } + + node = newnode->next; + if(node) { + oldnode = node->prev; + node->prev = newnode; + } else { + avltree->tail = newnode; + } + + node = newnode->parent; + if(node) { + if(node->left == oldnode) + node->left = newnode; + else + node->right = newnode; + } else { + oldnode = avltree->top; + avltree->top = newnode; + } + + return oldnode; +} + +/* + * 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. + */ +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; + #ifdef AVL_DEPTH + if(L_DEPTH(child) >= R_DEPTH(child)) { + #else + #ifdef AVL_COUNT + if(L_COUNT(child) >= R_COUNT(child)) { + #else + #error No balancing possible. + #endif + #endif + avlnode->left = child->right; + if(avlnode->left) + avlnode->left->parent = avlnode; + child->right = avlnode; + avlnode->parent = child; + *superparent = child; + child->parent = parent; + #ifdef AVL_COUNT + avlnode->count = CALC_COUNT(avlnode); + child->count = CALC_COUNT(child); + #endif + #ifdef AVL_DEPTH + avlnode->depth = CALC_DEPTH(avlnode); + child->depth = CALC_DEPTH(child); + #endif + } 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; + #ifdef AVL_COUNT + avlnode->count = CALC_COUNT(avlnode); + child->count = CALC_COUNT(child); + gchild->count = CALC_COUNT(gchild); + #endif + #ifdef AVL_DEPTH + avlnode->depth = CALC_DEPTH(avlnode); + child->depth = CALC_DEPTH(child); + gchild->depth = CALC_DEPTH(gchild); + #endif + } + break; + case 1: + child = avlnode->right; + #ifdef AVL_DEPTH + if(R_DEPTH(child) >= L_DEPTH(child)) { + #else + #ifdef AVL_COUNT + if(R_COUNT(child) >= L_COUNT(child)) { + #else + #error No balancing possible. + #endif + #endif + avlnode->right = child->left; + if(avlnode->right) + avlnode->right->parent = avlnode; + child->left = avlnode; + avlnode->parent = child; + *superparent = child; + child->parent = parent; + #ifdef AVL_COUNT + avlnode->count = CALC_COUNT(avlnode); + child->count = CALC_COUNT(child); + #endif + #ifdef AVL_DEPTH + avlnode->depth = CALC_DEPTH(avlnode); + child->depth = CALC_DEPTH(child); + #endif + } 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; + #ifdef AVL_COUNT + avlnode->count = CALC_COUNT(avlnode); + child->count = CALC_COUNT(child); + gchild->count = CALC_COUNT(gchild); + #endif + #ifdef AVL_DEPTH + avlnode->depth = CALC_DEPTH(avlnode); + child->depth = CALC_DEPTH(child); + gchild->depth = CALC_DEPTH(gchild); + #endif + } + break; + default: + #ifdef AVL_COUNT + avlnode->count = CALC_COUNT(avlnode); + #endif + #ifdef AVL_DEPTH + avlnode->depth = CALC_DEPTH(avlnode); + #endif + } + avlnode = parent; + } +} diff --git a/src/avl.h b/src/avl.h new file mode 100644 index 0000000..880e3e1 --- /dev/null +++ b/src/avl.h @@ -0,0 +1,186 @@ +/***************************************************************************** + + 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. + +*****************************************************************************/ + +#ifndef _AVL_H +#define _AVL_H + +/* We need either depths, counts or both (the latter being the default) */ +#if !defined(AVL_DEPTH) && !defined(AVL_COUNT) +#define AVL_DEPTH +#define AVL_COUNT +#endif + +/* 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 + */ +typedef int (*avl_compare_t)(const void *, const void *); + +/* User supplied function to delete an item when a node is free()d. + * If NULL, the item is not free()d. + */ +typedef void (*avl_freeitem_t)(void *); + +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; + void *item; +#ifdef AVL_COUNT + unsigned int count; +#endif +#ifdef AVL_DEPTH + unsigned char depth; +#endif +} avl_node_t; + +typedef struct avl_tree_t { + avl_node_t *head; + avl_node_t *tail; + avl_node_t *top; + avl_compare_t cmp; + avl_freeitem_t freeitem; +} avl_tree_t; + +/* 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) */ +extern avl_tree_t *avl_init_tree(avl_tree_t *avltree, avl_compare_t, avl_freeitem_t); + +/* Allocates and initializes a new tree for elements that will be + * ordered using the supplied strcmp()-like function. + * Returns NULL if memory could not be allocated. + * O(1) */ +extern avl_tree_t *avl_alloc_tree(avl_compare_t, avl_freeitem_t); + +/* Frees the entire tree efficiently. Nodes will be free()d. + * If the tree's freeitem is not NULL it will be invoked on every item. + * O(n) */ +extern void avl_free_tree(avl_tree_t *); + +/* Reinitializes the tree structure for reuse. Nothing is free()d. + * Compare and freeitem functions are left alone. + * O(1) */ +extern void avl_clear_tree(avl_tree_t *); + +/* Free()s all nodes in the tree but leaves the tree itself. + * If the tree's freeitem is not NULL it will be invoked on every item. + * O(n) */ +extern void avl_free_nodes(avl_tree_t *); + +/* Initializes memory for use as a node. Returns NULL if avlnode is NULL. + * O(1) */ +extern avl_node_t *avl_init_node(avl_node_t *avlnode, void *item); + +/* Insert an item into the tree and return the new node. + * Returns NULL and sets errno if memory for the new node could not be + * allocated or if the node is already in the tree (EEXIST). + * O(lg n) */ +extern avl_node_t *avl_insert(avl_tree_t *, void *item); + +/* Insert a node into the tree and return it. + * Returns NULL if the node is already in the tree. + * O(lg n) */ +extern avl_node_t *avl_insert_node(avl_tree_t *, avl_node_t *); + +/* 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) */ +extern avl_node_t *avl_insert_top(avl_tree_t *, avl_node_t *avlnode); + +/* Insert a node before another node. Returns the new node. + * If old is NULL, the item is appended to the tree. + * O(lg n) */ +extern avl_node_t *avl_insert_before(avl_tree_t *, avl_node_t *oldnode, avl_node_t *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) */ +extern avl_node_t *avl_insert_after(avl_tree_t *, avl_node_t *oldnode, avl_node_t *newnode); + +/* Deletes a node from the tree. Returns immediately if the node is NULL. + * The item will not be free()d regardless of the tree's freeitem handler. + * This function comes in handy if you need to update the search key. + * O(lg n) */ +extern void avl_unlink_node(avl_tree_t *, avl_node_t *); + +/* Deletes a node from the tree. Returns immediately if the node is NULL. + * If the tree's freeitem is not NULL, it is invoked on the item. + * If it is, returns the item. + * O(lg n) */ +extern void *avl_delete_node(avl_tree_t *, avl_node_t *); + +/* Searches for an item in the tree and deletes it if found. + * If the tree's freeitem is not NULL, it is invoked on the item. + * If it is, returns the item. + * O(lg n) */ +extern void *avl_delete(avl_tree_t *, const void *item); + +/* If exactly one node is moved in memory, this will fix the pointers + * in the tree that refer to it. It must be an exact shallow copy. + * Returns the pointer to the old position. + * O(1) */ +extern avl_node_t *avl_fixup_node(avl_tree_t *, avl_node_t *newnode); + +/* 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 smaller + * 0 if the returned node is equal or if the tree is empty + * 1 if the returned node is greater + * O(lg n) */ +extern int avl_search_closest(const avl_tree_t *, const void *item, avl_node_t **avlnode); + +/* Searches for the item in the tree and returns a matching node if found + * or NULL if not. + * O(lg n) */ +extern avl_node_t *avl_search(const avl_tree_t *, const void *item); + +#ifdef AVL_COUNT +/* Returns the number of nodes in the tree. + * O(1) */ +extern unsigned int avl_count(const avl_tree_t *); + +/* Searches a node by its rank in the list. Counting starts at 0. + * Returns NULL if the index exceeds the number of nodes in the tree. + * O(lg n) */ +extern avl_node_t *avl_at(const avl_tree_t *, unsigned int); + +/* Returns the rank of a node in the list. Counting starts at 0. + * O(lg n) */ +extern unsigned int avl_index(const avl_node_t *); +#endif + +#endif diff --git a/src/avl.o b/src/avl.o new file mode 100644 index 0000000000000000000000000000000000000000..baef9fb812426584742d84fe46a19ee0aea47f62 GIT binary patch literal 21276 zcmb_^34ByVw)ef=eUm0zCxoR1AuT}`C6H)FMM2X9!tJC3351vcA%@U_I3bBihs6aM z6MQtS4lYj}bsTgSm(TrEpSyuD23(PEMp67moab|7T4xZIVO+@jpQ=;0o9@7K{J#28 z_bheJsZ*y;)vbH`hPS@?{8+N4>2?hlR0rxFx_D4gPC>09`UItBj%+x@GKz*xjcPjU z)z=3D>w-)G)Yn(dsi}(63iD=4jO>v~vsBYY5X+%$B%n%t{ffZyu9fx8txavjtFNDZ zv3Egi9TN_rBxrF$AcD~fy z!&)C)BmHGwb!=%QZeB8_`NVwDsdRKhh|AS<~pc7;a_`0W4yfsarXB7KzqUU z?SqcjaiR3YMT->WE4G(EU(+%qE6%^hw$<<}Sl`vUwz+jheZ7CS6q8ps$7>sFInhc} zvkmrmsl?d}tg3J7z-&nU`1SjYSMs;<{N9hH$f#&Y)wESoLS8A(BDGRd1cN$rqnZ^G z?yA(ZDvT9+8d##x(5VhXB_Ct$gZx2dks10gGdw@Dpx4~AJ7^Dk^XvMWnwqA%a(aZm zXz%VvdnC=3Q>2;Q{hon0C%||EOzL2)dBJ=Sn!5f*Jv~dqrq1@Dvw`iA z-F1DgoQ>o#`Rl#dTkWk{Xm)>==?nb>V*cy9_W50N4x8cnOfyus&)l>>=<|nHx+BSE zXwE(})Mth}C;3AI-hn@C>RmT_!%xg`vB%uh+jSsP0Pl|Tr`Sym{l3uHX9iFrZeM>A z^&GK-g)lj(MflgfGjoY|sdt&Tp4^MZll>##;Lx0cBjnF)@&~x+E*k&i92>e}v>9IF zLBy>20S$llcfCv9|H&UZWcAo2=10sQ;+JFg&lULv=p=W|e!3lD{?@w9(8UNX=E{{L z<5#`V7dm2wzI}NTBt8{ln|jGa5t}Gt6VwC~HEw^X-Lu5ItgcV=H+WZ3`*MDKp<}-A zSQ`7}^6&WJZ_*!Y?KB_0w1l%ge4v;?B; zKOl!tjhiSB4P zNoL4SSBwSfHO@%!|MOsIRAVt1PuyUv{8~gi%%&r*o(C|4*cCHn{`x{+`@&t=vX1yd zM}x!t7qq&Ass8IrJdq({Jpy(6!sQ;Z|2C zkYHxG)LpP2d-F3&)>{Z3HmBX*X9X)+sE2k>-{R6uFEpZ!(Hy@DY45j7fIgmviyO*^6*vg z7b`R9@P?80ecn)AzCXM`3GOAq#Srv|7Yo7tB#5-@3oXtsNAGN~m*i?7=L`Eha|%AT zMn$(By$?NlL$&!qV~*?Trb5}F;dMPdx6nI%r|a4hj13~Lp0kjA-Yv@t4;BIG{#JM0 z@U93{rH2Bw*9z1u1gesP8V%IqLSN`@3RJ2LlpBGPap6E=A4ITZpiBg6Hv;9Bf%17M zP`)1#s7e)pm_XIW1*-N0fvS{&;!)Ydv;OBm{T;beQ>fFyk3o#W@S(rbtBQ7dCyrD&iidV^rs8m4CFA_aI zDmy4Dm4ioRO$i!ZJ(rhTI4 zjOcu&(*lHyGyxXuk7O6@j?6FUJ-8Vf1^eZGb*zsT+FtMGxbj-P+&U+?C{g)~l07_BHo04R@OeTTGu zmat9_1JraT#l=P@dm|YIy|MS7W4)V^1)os_dJldm!9D`}BS{5&BNr9y|G^QNCn7SB zA~IJ*q_|-3u^k+aMaTIY6OJ7e&fegBayT-d8cY$3+3>)odhgv_dt5ka0bsndD4Ifj zj8Gp({1Hg&2+q^cMl?4gtXUwE3wp)F~2`je= zJ93k=JT$={3hXh%bM_qHLEp(bZi=24W5)mTc`z*Q`@T@WxDUB{4uj?6BOo9Q)4fJn=j9B{ja7y83yraxSPTW_U@b|2iXN^y^Ze8`&_ z8AmrPaHr5+OP+dpvsMK16?Q}DJ~Rs(x#kbODVx;!jn33x`a%Qz^mV#lJ^&&`MF4!^ zxp1pDXdqYvd(=}19m6JzsjCBnm+@+=u;h*P@8zFQ_p|$d5)9eq|bACKhCv5 z);HZ02ZHlxCVXai9?x1Bj{^k}1CJbN4`-!Wh6mghBabm;Ov^1@l0E!P)u!vCaAEhW5^wNTgUQO>S^*a3y zH0gGI(LI8^k2Gx_T&Ofi0O9q(^ah>IZ@PoHLcD>*9l0PosiA>v8Fe82;OmVo3j;{PShd)fy)L$-p3_G$oIJ1f~uu#Aq&^)+d$}g>OF#@(WM>` z=x%hGDtx>4FF=##lH9+MUrE%VzL{QlCJFaKWJ62)B`$^|{UHz;kK^Kc0hb5BwD~d# zOtKYcJqbuQY_n~-T2>ka1U1K&>&(hhREaG&GmAQ*E^Kpcx$dlqpxz?NwB>rV^h}_x z(YVl2z~;>)aGtF&i-LCxiT#AkQ_ux<4w-*bQRkBRKPu`xGXK7!ro(*J0Z_jnE9cvC zCuRK?D2l&rhAlTg>nNyuh`Pv@Td1Xz4@1);U?w}#jsZAz*Z&6Z6!>e)OMeDvwutaB zil?|}>2CrZ?s^|ry7So1O8*e(2z?m!sijjulH1Whh31WlUR{nCa?WvRTSCYjAWc9< zUWSXC(txI!%4q**B_B2q<8JV6tdN^_+@L9?Z$Qy$o8wy3CA_+1trU zgq4Mbb71Jz+aT@SjK&9n40#F2Q-mbr(o4u0xV(nyChL!~`YWo>x1E=OEqc^OZ>BDP zZ*}>5>QXM7S7=U?#Gv}A>!h93x8Q20?qOp&$;;G#_u84`((>azW#}& zQ_({yg%Jw?@b`&y{3)%EgAg~@pw`o5bRGR5 zMRbl5nM01d$;}it2rq25oe|)b6PuV=5w_4 z>IgYC6q5x%^l^Nv8Y&c52T36=>o0P$7<2=bGjrx8G&6Bb3<@daU27`EbLFBE%OCBXK_ssOZ(1)Y`A> zxl_^GJ;dsl1FHlrZaDp7AQOV1Y17tUmGqaO*OT%tN&gY&)2oTzE$LJw>$OC`rqdm3 zYzdG16YhHAjWDp048Ef~9|V3k>X18Oc~rQ z$$1R2BG8*y{57SPuRFiSuPc<>B1M$I@07q{=SyhkP*0@o`v8pTXX(y`;ES927$9^r zAMy?Q5}WrLaKz^Q3m|gyPN7?@eo%~GE$xFe2dTo|L2*Z!L@(>~DV$CEWD)ZzBIc7t z;!F|xXNaWA6S#Jo`u#YMf8)}tr;-&fTgi1JsKSs}PToAwqED|*ot_~K&l1Ka(w>q@ zj?UI|C@XU4?sg_Ej3VMjF!UVS1#{>(u^if?CsN$^P-BVi?154_CIirZVoL!ZMRt|I>q-DmRU*tqFhOB{6*fhfGodEJd>=R>%%1}hVLls8GR$yLo089w zIAld+IvK0NVpc1(O2*;}D2ha`662rnIB1c`izNLA&^M88jikrGOFEC}7fZScG@VEk zeHUmtfhc+}Xp!jEBGI*E<k&#PZIX}r-=d4j)j5Z@$Q%+w-6Bk#A+*CnN;|EH-Fe?~ ztCpo1fEEUB)v13E9eJu)_fBQwHr+WLIGuvTxhFy3QRfJ>i*wIX zAYx9h1|sJ4ULbN#Jrtt5K~mYpP=|F-CP%FwOcAxe0ufQuPCLMY>M297UO_{;gO`O$ zggeB5#_3gJ8QdYKm44cw4A<|J^shjRH798;{S)Lx1|_+sXX&SBWE=L3%%ss7A6=c6 zg_{8mQd7Ynsu!i$GtN)FAa!DzFZ&{>*bU9bEMGQmV<3~H1d&GJab!}V7Hvyb1`8nQ zBt1fgkOmj&IsoBRqQXR`XX!~o|wL{m7I2*2}1>aoE3C~2A*pL zU6A9WUHG!8lhfm^Xxap8VuZ>>>UL2Q28d24u|9W47AAA#h42|Xi-I_lf*>>@HYL6z z8nLzuH*lLakuppFHX|KZ9r>i8>QnL1yQXDoBS?X)A}*4E(ohup2byl!Qpya?`6%E# z!?B~xu-Vs^K4REwjih~sZco{Y#$mG2`IuolWF*~aIBqm-&UuD*i(#ZVPd99j8b*&{ zTVfcijpTMCX`7Kg-$;%^VM3z0S865aI>xS ztF`uwQs>Ht4ec)wy0Okk=`!pM6a||-@D5tXmd^AWMz>)nP7MqQd5n9+F2gp{UuNiw zF{o8h<%LF4H$=-w`MlCHD6cf^#ChrJrl@v4w977(+V)@ScVoB<47)S;6+a{oU+H|# zmYZ!P+yC8;z}hB`0=jMnVsy$ihO^s9sx^!Y}S>BpL~|MqIz*}CCGC)9HpbtC(5*x%J`#kv#0v%gG2_N=Xg$k|J6 z`8v|%eU0YA-sK!>U+MJN8=dbsPqQ~+cSHg>?C)v7&|KvF*qIIXMY{b>%{C5=AkHLk zifx)Mn(}n}3a1A{ENMa(9Mh1-_7ATT6Yadn7;-aOKH5-&l>3$b4u<~mh8ZFYe}kk= zwcn|4Mf*cnJ0Eh6yV;)6p{;2v`CL$Y7 zTeM(XyVkh6UDy}jo^?uvS1aA-ag6FbU`qV~~x4Tpnm@>yQy_ZWWw{reL5$ zh%~MWH2$=Hc>^Z9F`&t?b_5h1V=o$G@?QaZgwo@_N5H%8kkWXiM2DMhv6V@sV4gUhzdlWkeJrba1bpet)N}2 zNbG$N>&ruiFs(QWDv>{1T~#S`x;UDUR=tu|P-%T`X;mte99@^Phi*;MGgCY3>hp4H=4JAZJ8HkD`MCZpD zRGDwCaT}2BBX0PQbGyl|ja(46U~eRVwrji?vYMo@yCN7?%=xJE|o+ zjb4>!%B?}hrQQ>lS{)_qg^3~wlQomIqgJXeX4c>Fa7<;8P$=nzR^N3&Mz=*-y?8f3 zWj8}zqu+gxE9+SLI!kMTLZe!h2a}~z65eT2(Y4CoF#1%Kb&%07qO7BgZor!pRtz)x zXq5FQM!$)&qSVkI0#e5&chhAUZQ zu@bfZ>P2>hYc)4TJ4$d*X-7p8Hiwl7RubG@%_5gE9CT*XZyR5-5RSE$r;MndD$ccr zU8T^QELOtY(T@XE{%VON?5|4X&Zvk&pSD=F3O$LTjV!H1b#1{eX6>^KEl{X}*lnbX z(Q6M&D{XA1e`Hp)d;Ce$@kK=;+I>Z$RziT38&Sgv+nD+#M&sWJklYIOSgeE-Y$2nle;riBu^x zYA7Ma2k9%`@DC8F;XYSdeJ4a`1vf=~Nm!AWv&aU9gD$s(pNJBV_)6s0mN(1pW7PU# zDRGf1k&}e@6Si{5GFQ-q@vec-m1$PKEP0Id%Sns^8zOTM_SH@#9hY>PD)@t@PvNm7Kb(HgY`F-8Eib*^iIYkpID-c=@sapXAMj zUXxyf&?Co4R89{Zb?CU0q|u*q0u?kZ4N6NotyA|yQHP)w2)vQ&cCO@lQpaf;t0rDQXMqm8e@$>DTSYP4cWIKFp(dT*0bFV>@#Z5+049R@N|A&5+6Llt{3`LwxCiI z{g^rmKEDTk5o`|S38#-N$hXgM`GUtSJ{ke4MKw;P1b~ zZW7Cjj|0)J_CWFd0q7p?lRng-K1ZRy`OJF{Ht|oJz;SE?PV_WXD&liYK=raD69%aL zDSpaR`Aj@2gXQzTB12}WQl3+x>85pI$H*O6HI((ghRb+B6#BN9oS>jUioW>GCO&`L4EQ*N(L&kDEfH|Z41Kl3W{xq7J|1Qk zIMKqfAg$wyhbgaw3GMCN&L`u9_C4HwJ788St@z`Mdc%C_3EB^^MxInl(5J9aDMCF> zqrZhZzEp7g7|Bj(Z|8P86c6U#9&Vq1f_4w@6({mfeR*1Pg8tR}y*{+-VkJvKbye@l zdQZ@<-b%fSwe@z5kuW~BSZMwSZRdMi9yAAQSI=_PKT3jq_3TiF#s>UZgOG~3sAKva z;0ye*y z-!)mzdPr3E;`}7}yCA=gZOjp?0OzqtA2Jpza|m^Yc}Y^FbxXWG2C8SZ2_6mU4?(*a5U zd+b9%^L`Bo)x1%FIJ@S3lP38;10+5^e4~M>0Hya}3EHVfY-wCm7~39L;bu&PgP{kzpajDGak2rZQZE za~H|CGVH@yi|FSW-plX~hCgN4$Z$Ht0){W)U_<&(GrWRfC&L%hJV9+k^UzPf6wqmhW9c26~mhtZejQ{hU*zFWmv^9 zgJBZG1IP#JZy&=jL;715;`Ze~mmN z_&LKA5Jc0-gy4%vR)TasBp6^=$M8Ic`3xsA%wbr_kd!gyRSX?;#w41d<0#V%9du8j zc7~2aOfz(R#xz66r%W?+e9SaM#{s4pIzC{Up<^G@3>|-9nxTWvh2#%I$LmZpbnIf9 zp<^f03?02pGjzPbG(*QTOfz&m$uvX9V@xx2Jj66Z$Nfw*bgTo7C0@nQ5oDU7qn&An zjuxgFI#w~w(6NGPhK@^^X6V?+`Ev=wB8DD@8p8v~1InXbhT9l!WOxa~B8DD@8p8wV z7idy1!)**VGQ5Og5kn6{0HIe59(d3BDzj%jsnAEw#?soaT5^6K1v5z`02qjCt9+SLy}Z!oRS z;m4kc1yH+p0 zixXc63KX9R3MfAs6KLE=f-yWLs=m6D7NXzW)vp#G;BhC>kMD$B%-46~3qL6p{RmG8 zML)M=()!*mmgW!cVp;NwyI79+tSyFt{~NGA(9zM>Av=!!woXXJeq0yBU^k +#include + +#include "macros.h" +#include "hv.h" + +// Interface to hypervolume algorithm by Fonseca et al. +// +// @param r_points [matrix] +// Matrix of points. +// @param ref.point [numeric}] +// Reference point. +// @return [numeric(1)] Dominated hypervolume. +SEXP computeDominatedHypervolumeC(SEXP r_points, SEXP r_ref_point) { + // unpack R data + EXTRACT_NUMERIC_MATRIX(r_points, c_points, dim, n_points); + EXTRACT_NUMERIC_VECTOR(r_ref_point, c_ref_point, len_ref); + + // allocate memory for Hypervolume value + SEXP r_hv = ALLOC_REAL_VECTOR(1); + double *c_hv = REAL(r_hv); + + // call Fonseca et al. algorithm + c_hv[0] = fpli_hv(c_points, dim, n_points, c_ref_point); + + // free memory and return + UNPROTECT(1); + return r_hv; +} diff --git a/src/computeHypervolume.o b/src/computeHypervolume.o new file mode 100644 index 0000000000000000000000000000000000000000..f51bdb169265adbe5ccc7f2d2d2b94c833b824d5 GIT binary patch literal 3440 zcmb7GYm5_B6u#4Tx@EgUrHgGcRoitBTYU-XE z832lsADYN7*sk_IDaePEZ-p5Pp)bJpBe+db%4)Gbt!Vm`fj&hU+$`rdxx>umiQuB3 z9TpwD1bmdL&+$0YMil#TajX&rtL+BG`bN1x<=L z=KyM3KC5Q(I4|E7iVjeo7m1GtJEARL#VEDqyM`A6+A`w%QF81i*_O{N=+mmzy51hz ze|Nr*AF69{^zap`)rzXL@1HuOOqp;R&IIrKog}^t@r@84nGDPsS4YHsj-e8A7J|v&1$_p7e-YrH zfO`O89vgw)-0(_bZ@rSuQSJ4yjQx|mSss_i?R zev!3*;of2UtgD`i!i>)_2bz11@U%N`*glgv^!)ij`_`cSYu5gqk6)F$x&zBkHO>1b zbKs1$Ro*6Vmv_j@_R|uKIp>HPV-Yrh0Br&~ccI}{x5kmX4q^NSjEf2v5cnHD&K+p` z`#_v%m;vAX$Vs$fn4DMv@g+B}qd5=ay9M)8P=7|Vn+Jaj>QyusBeD25h(FL=f}=q+ zg>Ix_9msJMy6=P_{Q$V=7zjbS3^H^UL_iAP0El4WLM~rl2+73&FbFH6y`n!F3*Hm` z{>4bQ(1W{^1b>gPG31Xv8h$LiGBVt|0q%-lWWJ!|8tw(wwaiBmXo_GqjRx1Uh0uxNr3zk}rLjaFnlk7{TtvP*l`=nB17SN@3 z77fcvnQGOrw3=aRu$XDBmQqV*$}&st)Inob(+f4V+_bO(Haavh!bNUqfSF3wfJIuM znCg^cu#)R68}(vEWo2!a<&=u5D>VkeOn1jiO_NMzW)5b8z(62_<8azyLV=CB>SV(; zN?HqaH$w8Rx~t*Z<_%v(j`s-tOnE)v4O&Bvw;Ms*TZ8&(5Z?|qf$eaV1vyFT5ki5r ztinz8%|=K|%9&9)LwAd9h_a-!>oL`uL6X{+Gfl&K!CUizF@?^PGo{1qIiE6=o@P zXjpU5g#`#dLeQ|@pbJY7xp3Ffg#VSP>|o+PK&WSlCFFm^!C;zjCE+r{eS~;{;rH@{fjro-?|H%iev8p2 z3}neB3}nbAgj9ImLENY=6x8TI1fj5Y&kYKwj*8T z(ULT~ve`^-l}+EYyW6I({jy~r`MUksv}Hk#A$BMWPN3kDE-`NzZNZziAqf!e@661# zua!XA?w_mqeD9nyGv~~i-<)|$)>CJHdv2DZI0_X-nTf}R=OBoFPASUAz;7*3lyG>Z zzi}m)+~CVekUVrMGw@JMAsjX$TMa=Ae!9Mit3+PUd=V4xoGMPtPCFY8$Bg*;*kn?! zea*Lttktpv@s7#Tic%!?QJX0Uec_HsS8pWB`RVpOB=!A0FRduV1D2_9o_(!7-La_# z=i1l%IU(wk8i=3P-YLzT3Wr-`MzphggG|rWcR=cEm-!ZNshMUg62swm_m6 zdk^QQx7U5Uko;SzfcRV;GK)vyaHJ#L9&One5vjTM9hCZXDMUQ^-1`OcyHRvlI9#i( z)NG}Ek|sx`;Nx%MxFE^jNH6&wg>bmLD;(RlaeYr$xYvljuRcp1>9gBK{-Kg08N%TW zJyX4UzI`sAC|D%R6Q8TESmsd_4tKSDkbQ2*Q`jmZg?PI@YVRd7kD_q6t;J}WO31}; zuTSJzzDE29_v`xD)F?Z@z4uFfR(nai-9G!6qde+Y@_}$T7+7+LzcygaDfN_#;NfpY zk^IU_1yAj?>$(`c2ah`^J<=L2i@~9Px1}d~ zPpo`-86Swml;0VN#yT+*m4CLZvV3_Y))SAmMq=e%o$KjYUe>k^5%iJ%4m_5AAORQo zhEl#_Tg-@TyiheYJ<-Ubo_IGL4mDh;POW@;8ClQ5Y4!X&>TxLL@mSQdN>o`^J!u@u ze-aPjGrC++N)a!Ryn_Cs@F_eTKkz##As_MxpSBalZR zk3b%QJOX(H@(AP+$Rm(PAdkTRBLo_B_Ns0!bycSH{v$@AsZ}%$wzs#Jsq_0d3OqkC zg#650b-qv06GNWtlM9~I&CAvKwfO9YJ*W3)jmvfOE_dY-z5f_=EUu_}GWyzBh{M}w;Wq;7^?I%^u$`;Qxept-@FnWZ!TsLlp; zGgca8SuOiae`@O$J8sm?8js$eijQR~Vcj%)YTT4D6kt~$%pzpmfg$=(P{|Y`2{08= z4!1AysB4Y3R$Hfq$-F5%rF~=@4Dwf=B717d9?-nHDg37ycFB$_baRadh)wnnc-V71 z{v1{RLXb^N*4WPbhxZ@wo39Vm@y|mL)8}7`K13g(q5rnHb#@n^MPI4&1%7v~2(VK+ zd*>1QWF`IOs|f29sqje^KBi~wOcKr0AMII`GSQK3(iY`J? zX+L5KA@?0|binih7+rt@ocYF4hw}+=_9PZf-8hT%?;2K5aBPys;=#^Fd-Csxj26tI-%w8lS(~8Z+r%4K%~_pHln&0X>LW z(Ixyh!2S?0VvDi0kZ+) z_!u$O-UWk+cd}~I0|(R1#qP@En41q;Wi>C-FgYF3e#Vn#adCE}YE(^JjgNu@rsfGS zKh08xGsrqzNWc)Uruocp>Z~5d9SgBePl$yuKQ41qsl~1k!^E}86*A-QUyx&m&=)f` zb6rJtFl26ry%?E>Q#0ZAgBok92%5{Sg2$-fY7`8bt9ilWR1j`A$W~W`P`gx^qH>KW z7cc`Jf8{Hat)gYeF^C#9wyMII;a4AQ_lgQIi^x z3!nlX0xIw)0JY2lzy{Q+9H3Uc4^YblP_k8pX7wKj>Mfwbfhx6tS^z?>k7$cKgq7Zj zhPi_5Zn#PVya2;DioS0WD{BaIGP>#xEUf{viTmnttn28}0gZ*w_0?86-laOq1>uwa z%41M1@zv&wN)>QOGp9X%J7 zLxjq*3!}1m5i-@jO{h@Yy}~>G00^L__U(ZM1d)jteyR4QC_(d%{~R&KrrM`*8MNqT zOQs71DQQPsc|3Dz<&n(N%GCHSNK_sd^VPXQ8feF~T}@sw8_|ec?WML7RpYH1+&l*b zj23P8YHy~OGFEH5*LX9L$|Kk>oEsqdl(sA671I$(>5xCdMqNpLcva zH};5ew+7QQ(aO|0i2A+bYlMZxFl&1M3;R;LR(q+P3Rhb11pwqn(0)j0{|wiz@`fKV zK1H|`kjhkMUS-O@{+vtgf*1T3A&?qBA)--=j%Q|89?RTbdHh2ovV+>2F^Gd)N?Pht$54ND9a;AYy5x*aK>`2@J6f zr2$g_w<2VAy9K_JXN%Xh?5_k~UxBJgyfT}n`}=u#olT$l15O!@L3UbCq}jn@$W zE*mCO+Ez~Xm32PA7F9Q~?b5ORJYryyEC!n?)Y;|akJyYYF3Mn|u@qanUfKy6bM>l4 zZeuRr2F(s8>s%Qp?Sr5JPrZk<+(4o?>!zT)3^^=*8Yw5)V~X`tUP^P z7Ze9TxM=!g*mefYI+&F*W&o`02%SRu0WFSK0ugzn!}-_Jr%m|}geOmR7(zY#RA-(4 z#Cc9S!v~({Zn%JHJI-Cmv^VgD)JfAWFurL~%o;ME07A<3MF*296xaH|J3& z%+^hf);!*QI9!-ZXbVK$r@0M4xyjRnlX4xXEiO=Ole=kf@APmU9yms(XTd-mv^kEK zg4Wm&L0Z7c$Aawj^sDn|upK7yHIR?sQ69zwuB4BG)Y(w_MINKl6X}B>1ME@V+&)sv z&egIvg6zN3pFrW0np*5Ed4$XT8!s=;AX6SZeZ1p_K*z*(zs7#0(e$)^IMn+3;^eKb zLYI>MmZ;I+|H_Ulv}9dYPhR#IjXgbi<7VEV6X|<-k*0lm^0ue(<4zpTq}YrQ z29cKZJDC4#dh2NLp5HBmBcLY;Yp;MUghPCgy#qi}>6dvDrrR%+a{g)XL9^*V z&}@LJ&~Si_2AHqH-~ZbkZLkBi9?=J0EDR(LFAOB>J8IdRI(usBycS4aS*5epBSX~_ zWGmR4!Hy!n1@TyrMbn^f2Yp&kUb#?bAJ-jEkup7TXkqrr_#=TKl2Zb$kHG_m_2jLO z5u3hMm-QW7|4S-@sZ38i8BR67?8RJc=W$GOcIGh@5XVn=ilF-QunOqz2eU@e=68I5lGfy?)4uK%MhG#u$HA0V};Cr zkUC&;zs@`KtI%c4No+r;90a8p3qd*rc`e8SLCy`Bv-E^=nqQU$%(+8)!f{&9rfj5^ zc+VNtYuWgM0o`mkfDw3zclT6%x*1LDiJ7MnV{Ya#8b>%9CrA3Qb?Te>5?{h6?_l}*01qHauFhi+(Kh5_)W<7LQ=;m^i zHC_TO`OPT5@-+grD0?J6dWY%6?Vb;XphPfBQ3BdS$+nU#CKdIsM=+byvU@#jpz(+S z+l1ZAV7;HcX@1E;K{LJ&!{mN>{^R}gSHR(9wd!&X__4F&)94#AUpK$#T5R4rq?_~GtJN~+;Ff0E zlTY?9rhlvpxUDk1FX*Oo6m7(E;6-euQyS%MhLbp7(6fVD)${S;?BiN?z+W|f?{_k) z-Cy9jfo#fObsBT9+P@C2nwdj0y0p4oeJj;}jBqiEwd_-v12N-=*hHGXZ+}uP#*AWq z!JST{XF(LDPYA-5epV2r>HiQ!Mfy=9Fo!NivBgDJY>p2TO``iJ!e^ZsCf{fEX*23f;?Y_4mo zK2nOR)#Q6#yrEz^36-Rqxq8ORf|I_F@7(u~3P>nmAjaerX9t zsh^DmS^OZ~Fu|8^3a}@G>J9pE2Sm%E5Mf9T={$Qn(Z&2Q5^!u%c0n1$M2`Ujg50k>pGte#0L5yP@ zz~cv~QfglIaWXgAM0+dt_QFd(yUaQ7Ws*R7GRGVIwI;6Z1dje|I(%;07 za9t!*Pclb&ZsYLy2`u6SWx&LZE}x|oWde)aOtN*){kTR6n74!GQiWMQ~OGw56Fj7x_R4E+{I&V@D0Q1&hcy(0u94B zTNDM&W?Tq*OY~$DKhKf-1#fAF9v29fHZUvWh{cCB+3cb!)(&fldrLe8_ZCxiNc5-1 zUx4Fr{EI^ecs;21E~=M3pWgBsCMzs%{iz-Eh98C>a0}Dg!Eo%C# zfaxrO?U)LfHkQJI!YP?iq*AAa*qdZKJf@fH zr|a#Y*hsw@tPu3yw*=|mr7&^ZCy1?O~Y&L09xc9rYewv-9}mAmJ*n? z1$~=7ds13ts|HmB*nh^qk1{K$%(GUR#MUCk&_e7Fn9vSFSHo~-BduVQ^93DrJ-Cb{ zsQoEXA{JyBL1`konao^tzRA!lt!^Nzk3x>>vFnF|nj+DVwLaNgIHC}1Jc?|8ihKb? z=l~je@mDb&Y^Hqlb0fzg#XVxuap)@H@Oy;8Z5wd)ot}-(Ch64p_s|Zq6t1kUodpON z6&<9ek;6yVg_6O?uy@8XpEM2^`LWm27Njb(goX_dk`JCa4F&v7T~Zu8cbbx&XOK#- zRb-a)9G2PwT)Y=(Rj1?VtGbHCLr(GqMZN2l%yn3h;AvYbAoC)oyHzElcsG_{VbUr< z6_;Qe{FJ}yF3Wjv&T3J1)krF`a=0#34~*!k-K|I|Y5 ztK*xk#aR~mp`ZQAhDV6~9`g;=aT#izMH5xEx^^xXzwG#q<@|QbRz6?`OE6K)#eOMJ z6l5j1zzx!H<_Q?!CBV*gC2rv7A--$?W!l`6`=B8=( zqcGKZLAYoFY1r>)umQpxV}W3Zy{`7teu^fQcQIic5Z%H(UjW93cjFUk=>Qh_g|SJ-%s%O zlXzniOkU*&p{>YXcO&Ax9h5kZ;hF8Af9D(j4N=z;b&ROb5;aUz4N(shMgKiF{!^kN zME#hkZA9%OY8O%8AnKb$eTk?a5cQWt{e-BkL_JJYH&IUz)kf5-M6Cve*cJdW54-`O z&=EXf7P5s2{CBVjP=N53h!9^FyZO4dn?3RFXr!gJqh)xwz3`$)E%8Ys_=`BzSJSeoUP;pkwP#hIF=}*j43Vo{L7{hHty-RU;@QpSBalZRkHDXdK=C#1W?VNYm-@o=6L@WrXtY;Zw#l{EhhK(j<3F+& zUXS02!;j+aOv_uCtYfdT$zALVFZXLf1+<6g_HavAS5NDmkyfK8N~xtHwJXxS!RVlb z3X#y=+S3&yexc-}JzF^M6MS#9$AA(lp^KDww=J!no4Wq0I~~MEHpU`GxVObHBGK;3Kc|dcl(8v>%Oh|i zIKuc@FlWo=uCP+q+Ox4YZbWK(Hg82m!QodWL%xYD(wJ2Y8bapF0cElnr(bf*-QAfM-o+IAf)!BVdxVxt< zqWr#~Ez%V+BFY(D7vZP23}uf4S6b1qLG3t=-~XDp1s6+(@*@XfFOq+U&jRiHS5f^xiI zeWV?~9HxvFw6p_hB zHL?b;)LOj8niQw9_b>SU*u*WK`sRI(-U^?mxxS))#F280T)ea3Ktak;kK}!h19K}T z;P_`>n2Qh6l%06!%FYU;qXbU@o<@|RbT1xCr~lxfbT_)ctfd!sjAiHj%vqt7p{ql? z%pCmMn)UN$$WjVPQ5#V*(C?br`Eu~>g%0Jl6)PI32q`PW%lf^{HQ?!c)fP|rWVi

)&5V-GP1#FM=1YR)R{h)CYzDWBx)MfRk;ohSKbDc>#SE#8uE&XKqCB%hMk{^-fza`}@-jc7#k+<_CuXsd({fMZ3 zizkxs^WmXD{1Uu;#2^q)@;(o7_^nkuByaH`KUDuUdr%Nz7!Pm1hd801>Tf0nf%H-R z&BP!OPxYrf#Nqe!@Q}R4gZxna@4AMh73Fn2ByaJSyd5v9-_G;)OZih$-Y0oWzMhmL z5Kr=Up5*sl!?Tr{ps9X~2g&Up5!sXScAn&$uN4J8Cgm;Ol27Hx+j)}rd_olXl$5u4 zOMZlA2?XNFemhU{6Ex5f=)3!5zr};(_H#kd#FM<8=knJh4dD*tae2vG^8Cxx$hYL} zJjr`L&Bc|wq`bvj@+px$DR1XVeoV@DNO_C5p&CEt=GZ|B>nkPNx=@7IynW@RG&F?gDP58OF@{`HlKLa!p9@V9u7y#Kx{ zbr4VOxAPxVKP?FM`isC{X#O9Q^?w%mRKLYr^;hKhpPeWCTsH~>ZEKPjJMe}#=F zd7qSDBIPaKlAplWA}A0~?YHwJUm-s?S}Wx(-jb(vkAfv{@pq&AC!P4*44V}SUjXmn zlz5qxkLB9>U0Ra2~+YA0%(BukT6TC-X1Hu7k>3>$h9-)_5q9yfwaVl)UBde#tjWeGQUNNxsd>m+Nhh zHz9$wdr%&Hwoj4SBMq=d~I*oGtMss;A z+FGvQ6on7lKU`rd+oOc(;MK}cLY!~PK^BI5fkb;KPAq!@q>v}P6w?p9dzuxJo4C+#NmghL?i24x>~whBhzGZ zm5NOJ>5vPH1EM$^*%G9&QwKG9R<#{bKQ#0g5JgXX;Nh2li{is|EQ&&#byntTTcQ!( zL*bTKtAZ1BI2vh<<9Dn#10kxG@&!C)(@*URE-a>dWb;^gptp6Et?$~W`imF%qs}+!_2>EFG$-%u=I4K|Q_x|X|H+gKIBBjAPC2%H zL|tgaP4m+abd{;D&+R*Yy_9>gM7+FGB2ru?af)9f5oKK#5j}Nvg!?UVX#_%C7ujyS c#GTRPZIQSqvTlgZzZ#NvLg!tP@%dBvA4(6pO8@`> literal 0 HcmV?d00001 diff --git a/src/hv.c b/src/hv.c new file mode 100644 index 0000000..4401b4d --- /dev/null +++ b/src/hv.c @@ -0,0 +1,525 @@ +/************************************************************************* + + hypervolume computation + + --------------------------------------------------------------------- + + Copyright (c) 2010 + Carlos M. Fonseca + Manuel Lopez-Ibanez + Luis Paquete + + This program is free software (software libre); you can redistribute + it and/or modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program 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 + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, you can obtain a copy of the GNU + General Public License at: + http://www.gnu.org/copyleft/gpl.html + or by writing to: + Free Software Foundation, Inc., 59 Temple Place, + Suite 330, Boston, MA 02111-1307 USA + + ---------------------------------------------------------------------- + + Relevant literature: + + [1] C. M. Fonseca, L. Paquete, and M. Lopez-Ibanez. An + improved dimension-sweep algorithm for the hypervolume + indicator. In IEEE Congress on Evolutionary Computation, + pages 1157-1163, Vancouver, Canada, July 2006. + + [2] Nicola Beume, Carlos M. Fonseca, Manuel López-Ibáñez, Luís + Paquete, and J. Vahrenhold. On the complexity of computing the + hypervolume indicator. IEEE Transactions on Evolutionary + Computation, 13(5):1075-1082, 2009. + +*************************************************************************/ + +#include "hv.h" +#include "avl.h" +#include +#include +#include +#include +#include +#include + +#if !defined(VARIANT) || VARIANT < 1 || VARIANT > 4 +#error VARIANT must be either 1, 2, 3 or 4, e.g., 'make VARIANT=4' +#endif + +#if __GNUC__ >= 3 +# define __hv_unused __attribute__ ((unused)) +#else +# define __hv_unused /* no 'unused' attribute available */ +#endif + +#if VARIANT < 3 +# define __variant3_only __hv_unused +#else +# define __variant3_only +#endif + +#if VARIANT < 2 +# define __variant2_only __hv_unused +#else +# define __variant2_only +#endif + +typedef struct dlnode { + double *x; /* The data vector */ + struct dlnode **next; /* Next-node vector */ + struct dlnode **prev; /* Previous-node vector */ + struct avl_node_t * tnode; + int ignore; +#if VARIANT >= 2 + double *area; /* Area */ +#endif +#if VARIANT >= 3 + double *vol; /* Volume */ +#endif +} dlnode_t; + +static avl_tree_t *tree; +#if VARIANT < 4 +int stop_dimension = 1; /* default: stop on dimension 2 */ +#else +int stop_dimension = 2; /* default: stop on dimension 3 */ +#endif + +static int compare_node( const void *p1, const void* p2) +{ + const double x1 = *((*(const dlnode_t **)p1)->x); + const double x2 = *((*(const dlnode_t **)p2)->x); + + return (x1 < x2) ? -1 : (x1 > x2) ? 1 : 0; +} + +static int compare_tree_asc(const void *p1, const void *p2) +{ + const double x1 = *((const double *)p1 + 1); + const double x2 = *((const double *)p2 + 1); + + return (x1 > x2) ? -1 : (x1 < x2) ? 1 : 0; +} + + +/* + * Setup circular double-linked list in each dimension + */ + +static dlnode_t * +setup_cdllist(double *data, int d, int n) +{ + dlnode_t *head; + dlnode_t **scratch; + int i, j; + + head = malloc ((n+1) * sizeof(dlnode_t)); + + head->x = data; + head->ignore = 0; /* should never get used */ + head->next = malloc( d * (n+1) * sizeof(dlnode_t*)); + head->prev = malloc( d * (n+1) * sizeof(dlnode_t*)); + head->tnode = malloc ((n+1) * sizeof(avl_node_t)); + +#if VARIANT >= 2 + head->area = malloc(d * (n+1) * sizeof(double)); +#endif +#if VARIANT >= 3 + head->vol = malloc(d * (n+1) * sizeof(double)); +#endif + + for (i = 1; i <= n; i++) { + head[i].x = head[i-1].x + d ;/* this will be fixed a few lines below... */ + head[i].ignore = 0; + head[i].next = head[i-1].next + d; + head[i].prev = head[i-1].prev + d; + head[i].tnode = head[i-1].tnode + 1; +#if VARIANT >= 2 + head[i].area = head[i-1].area + d; +#endif +#if VARIANT >= 3 + head[i].vol = head[i-1].vol + d; +#endif + } + head->x = NULL; /* head contains no data */ + + scratch = malloc(n * sizeof(dlnode_t*)); + for (i = 0; i < n; i++) + scratch[i] = head + i + 1; + + for (j = d-1; j >= 0; j--) { + for (i = 0; i < n; i++) + scratch[i]->x--; + qsort(scratch, n, sizeof(dlnode_t*), compare_node); + head->next[j] = scratch[0]; + scratch[0]->prev[j] = head; + for (i = 1; i < n; i++) { + scratch[i-1]->next[j] = scratch[i]; + scratch[i]->prev[j] = scratch[i-1]; + } + scratch[n-1]->next[j] = head; + head->prev[j] = scratch[n-1]; + } + + free(scratch); + + for (i = 1; i <= n; i++) + avl_init_node(head[i].tnode, head[i].x); + +#if VARIANT >= 2 + for (i = 0; i < d; i++) + head->area[i] = 0; +#endif + + return head; +} + +static void free_cdllist(dlnode_t * head) +{ + free(head->tnode); /* Frees _all_ nodes. */ + free(head->next); + free(head->prev); +#if VARIANT >= 2 + free(head->area); +#endif +#if VARIANT >= 3 + free(head->vol); +#endif + free(head); +} + +static void delete (dlnode_t *nodep, int dim, double * bound __variant3_only) +{ + int i; + + for (i = 0; i < dim; i++) { + nodep->prev[i]->next[i] = nodep->next[i]; + nodep->next[i]->prev[i] = nodep->prev[i]; +#if VARIANT >= 3 + if (bound[i] > nodep->x[i]) + bound[i] = nodep->x[i]; +#endif + } +} + +static void reinsert (dlnode_t *nodep, int dim, double * bound __variant3_only) +{ + int i; + + for (i = 0; i < dim; i++) { + nodep->prev[i]->next[i] = nodep; + nodep->next[i]->prev[i] = nodep; +#if VARIANT >= 3 + if (bound[i] > nodep->x[i]) + bound[i] = nodep->x[i]; +#endif + } +} + +static double +hv_recursive(dlnode_t *list, int dim, int c, const double * ref, + double * bound) +{ + /* ------------------------------------------------------ + General case for dimensions higher than stop_dimension + ------------------------------------------------------ */ + if ( dim > stop_dimension ) { + dlnode_t *p0 = list; + dlnode_t *p1 = list->prev[dim]; + double hyperv = 0; +#if VARIANT == 1 + double hypera; +#endif +#if VARIANT >= 2 + dlnode_t *pp; + for (pp = p1; pp->x; pp = pp->prev[dim]) { + if (pp->ignore < dim) + pp->ignore = 0; + } +#endif + while (c > 1 +#if VARIANT >= 3 + /* We delete all points x[dim] > bound[dim]. In case of + repeated coordinates, we also delete all points + x[dim] == bound[dim] except one. */ + && (p1->x[dim] > bound[dim] + || p1->prev[dim]->x[dim] >= bound[dim]) +#endif + ) { + p0 = p1; + delete(p0, dim, bound); + p1 = p0->prev[dim]; + c--; + } + +#if VARIANT == 1 + hypera = hv_recursive(list, dim-1, c, ref, bound); +#elif VARIANT >= 3 + if (c > 1) { + hyperv = p1->prev[dim]->vol[dim] + p1->prev[dim]->area[dim] + * (p1->x[dim] - p1->prev[dim]->x[dim]); + p1->vol[dim] = hyperv; + } else { + int i; + p1->area[0] = 1; + for (i = 1; i <= dim; i++) + p1->area[i] = p1->area[i-1] * (ref[i-1] - p1->x[i-1]); + p1->vol[dim] = 0; + } +#endif +#if VARIANT >= 2 + if (p1->ignore >= dim) { + p1->area[dim] = p1->prev[dim]->area[dim]; + } else { + p1->area[dim] = hv_recursive(list, dim-1, c, ref, bound); + if (p1->area[dim] <= p1->prev[dim]->area[dim]) + p1->ignore = dim; + } +#endif + + while (p0->x != NULL) { + +#if VARIANT == 1 + hyperv += hypera * (p0->x[dim] - p1->x[dim]); +#elif VARIANT >= 2 + hyperv += p1->area[dim] * (p0->x[dim] - p1->x[dim]); +#endif +#if VARIANT >= 3 + bound[dim] = p0->x[dim]; +#endif + reinsert(p0, dim, bound); + c++; + p1 = p0; + p0 = p0->next[dim]; +#if VARIANT >= 3 + p1->vol[dim] = hyperv; +#endif +#if VARIANT == 1 + hypera = hv_recursive(list, dim-1, c, ref, NULL); +#elif VARIANT >= 2 + if (p1->ignore >= dim) { + p1->area[dim] = p1->prev[dim]->area[dim]; + } else { + p1->area[dim] = hv_recursive(list, dim-1, c, ref, bound); + if (p1->area[dim] <= p1->prev[dim]->area[dim]) + p1->ignore = dim; + } +#endif + } +#if VARIANT == 1 + hyperv += hypera * (ref[dim] - p1->x[dim]); +#elif VARIANT >= 2 + hyperv += p1->area[dim] * (ref[dim] - p1->x[dim]); +#endif + return hyperv; + } + + /* --------------------------- + special case of dimension 3 + --------------------------- */ + else if (dim == 2) { + double hyperv; + double hypera; + double height; + dlnode_t *pp = list->next[2]; + + hypera = (ref[0] - pp->x[0]) * (ref[1] - pp->x[1]); + + height = (c == 1) + ? ref[2] - pp->x[2] + : pp->next[2]->x[2] - pp->x[2]; + + hyperv = hypera * height; + + if (pp->next[2]->x == NULL) + return hyperv; + + avl_insert_top(tree, pp->tnode); + + pp = pp->next[2]; + do { + height = (pp == list->prev[2]) + ? ref[2] - pp->x[2] + : pp->next[2]->x[2] - pp->x[2]; +#if VARIANT >= 2 + if (pp->ignore >= 2) + hyperv += hypera * height; + else { +#endif + const double * prv_ip, * nxt_ip; + avl_node_t *tnode; + + if (avl_search_closest(tree, pp->x, &tnode) <= 0) { + nxt_ip = (double *)(tnode->item); + tnode = tnode->prev; + } else { + nxt_ip = (tnode->next != NULL) + ? (double *)(tnode->next->item) + : ref; + } + + if (nxt_ip[0] > pp->x[0]) { + + avl_insert_after(tree, tnode, pp->tnode); + + if (tnode != NULL) { + prv_ip = (double *)(tnode->item); + + if (prv_ip[0] > pp->x[0]) { + const double * cur_ip; + + tnode = pp->tnode->prev; + /* cur_ip = point dominated by pp with + highest [0]-coordinate */ + cur_ip = (double *)(tnode->item); + while (tnode->prev) { + prv_ip = (double *)(tnode->prev->item); + hypera -= (prv_ip[1] - cur_ip[1]) * (nxt_ip[0] - cur_ip[0]); + if (prv_ip[0] < pp->x[0]) + break; /* prv is not dominated by pp */ + cur_ip = prv_ip; + avl_unlink_node(tree,tnode); + tnode = tnode->prev; + } + + avl_unlink_node(tree,tnode); + + if (!tnode->prev) { + hypera -= (ref[1] - cur_ip[1])*(nxt_ip[0] - cur_ip[0]); + prv_ip = ref; + } + } + } else + prv_ip = ref; + + hypera += (prv_ip[1] - pp->x[1])*(nxt_ip[0] - pp->x[0]); + } + else + pp->ignore = 2; + + if (height > 0) + hyperv += hypera * height; + +#if VARIANT >= 2 + } +#endif + pp = pp->next[2]; + } while (pp->x != NULL); + + avl_clear_tree(tree); + return hyperv; + } + + /* special case of dimension 2 */ + else if (dim == 1) { + const dlnode_t *p1 = list->next[1]; + double hypera = p1->x[0]; + double hyperv = 0; + const dlnode_t *p0; + + while ((p0 = p1->next[1])->x) { + hyperv += (ref[0] - hypera) * (p0->x[1] - p1->x[1]); + if (p0->x[0] < hypera) + hypera = p0->x[0]; + p1 = p0; + } + hyperv += (ref[0] - hypera) * (ref[1] - p1->x[1]); + return hyperv; + } + + /* special case of dimension 1 */ + else if (dim == 0) { + return (ref[0] - list->next[0]->x[0]); + } + else { + fprintf(stderr, "%s:%d: unreachable condition! \n" + "This is a bug, please report it to " + "manuel.lopez-ibanez@ulb.ac.be\n", __FILE__, __LINE__); + exit(EXIT_FAILURE); + } +} + +/* + Removes the point from the circular double-linked list, but it doesn't remove the data. +*/ +static void +filter_delete_node(dlnode_t *node, int d) +{ + int i; + + /* The memory allocated for the deleted node is lost (leaked) + until the end of the program, but this should not be a problem. */ + for (i = 0; i < d; i++) { + node->next[i]->prev[i] = node->prev[i]; + node->prev[i]->next[i] = node->next[i]; + } +} + +/* + Filters those points that do not strictly dominate the reference + point. This is needed to assure that the points left are only those + which are needed to calculate the hypervolume. +*/ +static int +filter(dlnode_t *list, int d, int n, const double *ref) +{ + int i, j; + + /* fprintf (stderr, "%d points initially\n", n); */ + for (i = 0; i < d; i++) { + dlnode_t *aux = list->prev[i]; + int np = n; + for (j = 0; j < np; j++) { + if (aux->x[i] < ref[i]) + break; + filter_delete_node (aux, d); + aux = aux->prev[i]; + n--; + } + } + /* fprintf (stderr, "%d points remain\n", n); */ + return n; +} + +double fpli_hv(double *data, int d, int n, const double *ref) +{ + dlnode_t *list; + double hyperv; + double * bound = NULL; + +#if VARIANT >= 3 + int i; + + bound = malloc (d * sizeof(double)); + for (i = 0; i < d; i++) bound[i] = -DBL_MAX; +#endif + + tree = avl_alloc_tree ((avl_compare_t) compare_tree_asc, + (avl_freeitem_t) NULL); + + list = setup_cdllist(data, d, n); + + n = filter(list, d, n, ref); + if (n == 0) { + /* Returning here would leak memory. */ + hyperv = 0.0; + } else { + hyperv = hv_recursive(list, d-1, n, ref, bound); + } + /* Clean up. */ + free_cdllist (list); + free (tree); /* The nodes are freed by free_cdllist (). */ + free (bound); + + return hyperv; +} diff --git a/src/hv.h b/src/hv.h new file mode 100644 index 0000000..515cfcb --- /dev/null +++ b/src/hv.h @@ -0,0 +1,48 @@ +/************************************************************************* + + hv.h + + --------------------------------------------------------------------- + + Copyright (c) 2005, 2006 + Carlos M. Fonseca + Manuel Lopez-Ibanez + Luis Paquete + + This program is free software (software libre); you can redistribute + it and/or modify it under the terms of the GNU General Public License + as published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + This program 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 + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, you can obtain a copy of the GNU + General Public License at: + http://www.gnu.org/copyleft/gpl.html + or by writing to: + Free Software Foundation, Inc., 59 Temple Place, + Suite 330, Boston, MA 02111-1307 USA + + ---------------------------------------------------------------------- + + +*************************************************************************/ +#ifndef HV_H_ +#define HV_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +extern int stop_dimension; +double fpli_hv(double *data, int d, int n, const double *ref); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/hv.o b/src/hv.o new file mode 100644 index 0000000000000000000000000000000000000000..188daf53595c6b6b08b3fb2f47989623ecdfe6ad GIT binary patch literal 17320 zcmb_k4SZD9mA`LhUU>O@lMfIIkpTjP)`Uc>yW19bS(kpe&9-irYHJ5rAfVQOC=^=7Z&)fKY%N$d`#<;HI};}Hr~Uo* z{&@F%+;h)4_uO;uyYt=)FMsmE=ffCN>Dw?^z&sW>h+!kNJiTYZFJ{iAK#95Agpk9gleATV3O@3c<^%_c-LF;3)(AO#k zwn)rd2iAvJI#;X10U`C?!175Xj3;ot}(9%=8-!RL1jqX|diI#4{4gKN2 z!#umYEBn3Ow<^U^h(3Jh(2YkJ_Y&(5`TlLx9jz$KTyl@ zxa3iyUl|~ITuL6hO500E;{?C-X#_9_tm~sd{bMt9joSb(03!SDz0BP< zi;eIemvR3A>Ji@iY3O%Gq-fI#86%u@8SsLDaFFQhIiRTt?{$Uud>T6TtOi3~ zvmbFsjMK?AD-YoZ#VvBXHs74{=4J9LQtF& zCuPM^ZzQw>%j9u=>}DD;pPz37xp=GYbnGcBK14%D{)l!XvdK9=GV7ob8C!p=Ug+4j zZe{VQ0qe_%U?+;Ka~hFpuNx7^2?~^9P6L#D;!#^WlEm@4(ZA0<{bZ=O|FFA%PucXd z8-8$3AGls&+0!2_n|>L4u->%-sd_FIJn414{$aH1D8(h1>F$3Sd*B@N^K$b7AK*nj z6FW-ezL8dRBKK($aY`g#B1IM=_v{GYuQ^ZdqZf52w1xlMxZAzlUFBZk_Eqhx-$wTf zw9L~`E!1|yzy5l?aL>9F)V=yYMD<@rAr{Za#QO{O{e2$uCAWFeYgX>^n4iy?&VOm> z-Jf9+%4<)0B6FP!BePC;A~zXb9|ymV1$I?@rY~}J4ik3@o|P2za0v^&8~((mKeP~w zr_Aj0nxWn2lJ&7zM>l!Q-+9d|G{1ST^I1LxUr({KFUZv_%6lR1>D?{5Fq{G$8n(cL~a4jw=8|D4#pf4 z8$8^%1BWAOB%iUfDMY~Caxo%=6^+Zavm0Rm!?o#M4*;`09%l50NjA8dh#&;}{9LvX>`0I#l)= zk=rzHWTn#+nU-gCSM||;M(O9Bxe|JuK+H2i?8qtBv|IG3pZb#ZRPg+Y2H@B|n#tb8pWOPM?IbLj0=Q6z!+U00#z&+4jZhlBZ zALfZT@{l`Bf8YZofl%a$ROHE(jU5iRq^=Ii78Der!p2E)tCsCMdvmSe2TdTD5=mhZuc%6iEJI5=$-+Sla?x5VIR&yAab z36J^5&_55733q5L=q@+^j16r!bXE4ATXWU@(tU6a0*eul>8K^q<~7e*nA-v~$j!Nf zU4~xabtQRy1pUa)!2JUXZquxREH4$qMIu_=mL5TKS(I`CNO=zin*Pc}#@eYJ^SqSD zA<8o%mN;D>>s<#v!;xaZ9!kBVXP?A4$Wx@UzI+%$I9JUS5%EXQxkv*Edm^{z5A5I|Z!#Uf_wd!BZ=J^vTl4$P@||Vos{`Uu zZvHR!8?uoDtHWrcx>aAE#=Ga7+PD~J;;Z84iFos{a?@~L@@QT&Zx1oBoH+vvJP)yR zj@jplWB}&uf!infj-6xe*!x4zTk&0bx_>X9cj{lXFfxr>%=>-W^ol%J=Q18L3Rfhw z6A{A!1HVZSvNz8B^)kikuhcgdtT}l@M5}EYB86WiyUyJ-V7$~1(emaE87|sDDtDEc z@By89jI9WkoA2vgbe^J3iK(clP?)?lQ$yRI%(S}k-zM^ z`7Qcazr~~!Tfe9XwGN9s+*jZ44LUE<0)`Hp?g@*X6ngndVUIj1Y}Djq9*s1zGQH$bBw##Qm_D*nAhvR;LYHGYY*)E^|sPw=lX)US;aM?V=2;|#t z(;fwk5^Iia0@Kb)s?0V))6Pq3u5E%-y9lZZuGr?;CXCb^z_t=)*e1A`x*v^6eW?B( zm(4wlK#6Ud_9UQ(#J)-DpOVy#q@I4l@wS^t{m&&;LhAoZQq!Sc`wb{M(b#6#Cgf{- zLETQ&EZc-4?Nv}qh`Pl#VH#7(hh%#%>R}H1DF97LYC!YHKy3x;8lWjWa;Y1E@`|Zp zplL~upiakSyIy^aqp4(*WsE_sbSQrzpfBS@0Q;A|f}^BNUl~K4^!rHJKwKxqR#{K| zoOQ&#??H5P(xVznfT^DXPZ|d9*NM4o*bM+BO1U-zkZdRmDQr(d&qktdr?62FIkkww zb{ixqh|0GpY=rUAEvB$70>!mdP}nLZbtl5c)OAp3p9P*nQ^UYjaziC&Gn-Kc=SWyS=wqxf2^5O5+6(r{kdf8L9t>YtQEK zx?qj+81XN1BV+yvI2Gk{CHD;QAK}V<0@X8w+yFmlH@S8Tki$UOm~LPzsROg+ehmEI zaOLiU{a+B$3FNmx#-9dn@%bJ0?8?Lak-&Yj;>fuVE=sJyG>?f^dr zSMC{9E?i^Tb*QZ}emc+>g>;gPq-J9iiUd(9nM-~Sd&)gL9(hWo1AZv?&}`=^6_i#M zF|i1dN+Gg|rhDQYAk&FlLnOOqA1OJ6sifg6hiEGGK8}{5S?94lMDt~>u@{KTt~p5( z$7QeM@m|N{y^hCw{hs%luu{1077zTDC2y0S_xKb}YKLch8s8&zlfI`i#vOT&I!}5E zeVydDMXts3N1hz-QEp6X#8%E-eV>fFHiPG)&J(BM+2knd_5mrK_yWi=ME)-)bMU;A zXVHh0q7$1z@-UvK*2%f_j1y(Ofa4r+NitgnQ*kJ{Q!t)%YTF3~BMCVKWCD<}tPnM2 z1QfbjQGHF}X`IQY4MzSkjoLBi{1FJ!W(GPO70mI1sa&Vzw!p|j((MJZjO(uBx|yq# zhD*Vj`~%_Dc*KLs=s(FJ{YY5{-e?`qaQTRT!h=3~4(PEIkRMz0deA&xKjrdJ*)3G9 zLCmSLcvI((dDEForZ*~hpNp~EkZT=G++DTp*)0S1q}kflH1P6CXoAzU1oolKU@^e98Iwl2bQeG{fm8 z#FR{BVn()V&u}J<&UkZETBec(wrv`-XLz!2Nt=;6Gu4I(wC90wb80F^fnE zHtuF40PGpW#Vt5nv1g>FT7;I)brDelP093_K4rM#$)1c^NmY68v`Jv54_AnlQHahY zJvq&j?SeL4g#u!b9Vc0rM5Yi`rD&EARY^2kh;kw)6~h`niCY^H??(Esi5nzxWFcab zrIe;{J6!cGaBGyxY)KZ2lz(j+8uQQNRn671)O4-b9`4lU ztEzU_(`r`8HaT0(zNBVp6V%Kg{M-r6xbB zrrA&1v{`DheMhHlmcmI}lVZ^~Tl71en#J2Bo?=@L{mJ%}XVnb*I@|Sd%XVG1nxV~9 z(`Tz0_p9lTsu{MKnQF2&UYmMJ&3H~tzXWhWO@81I90~Qpk&M^S{Jxs6Z; zv^{OpwyBwqLV4~s`^wkU%r2l!T2KqY=YX1SAGJlzouX!t)7NT+Q060aT)pJy4H>er`!Q^C0oK!A1s7xeGDN_w!gu&!f@%_g=+c_(R+cK?g1F5 zCNH$FpQWZhttPiUh<>kc)Lz3xWLsl-O0j)iu|KxSHbL2{rvC|AI_4i%({BZ`KB$dQ z)4gi)@3rr!>08xI+YKYa>Tqr6yPk!wEO=QfQb%byhtyoJ*7c%ifjaKc0(GKoLbf_i zD}GMRwcqUQ&fZR7xNXKLTfsFyRMTHVbPubU_K!9q3R$-NOn4gt z@{LS2{jj}5%hKBIDZOg$&tSPAn-6|2CgvGvy5?C7q*?m`(ka>g^G%q9Y=kC9hUP`K z<-wx6mUz|7)#|w4z$CH{I&@%pBhX`2uib^@QtX=*?Wa&V+BRi4X5e2s!Nhd@@I6}= z(%((5m0e!%l5$tOKhW0L((Jmiu(+_um0#P0Zxae`EV{n1sBqf!DX#pQ>NbC!3j$pF zhQb>Qi`v?or&z0IYGoU`Praw6rLAphz~9=^));IFH2T}7u4)WU_16ZbwgqZgTd<|o zSJ$}4-%OTRV{?!N1Aad(I3FRtAZzef*KxYJrH<%k`VOX*J~C;o4*1cb$zLC2fyPx0 zL55Er*09=^5V#>|T-6N!NN1hDHP{em)izWISg^XWiNQRpUDL|-^B2MaGEoo9Tot(? z3Gmze!BDHOwyvqE5wG&=nuJFd!a7)8OQ@#F&w|8fu-+2DNCN(9*51;@`1cNV41;8C zwSnqjZ3DZ9)i*ZbD-Fi=1^o4_I@H0MTLl}R=rs9*c*jrbctyU%?ANz8HToLbFz?89&3Bb2=8H>HKY=Go-!r+paVK61}h+xuvPU3P*T%{!BSV&Of@u5It zd=M5Uhq5k~=Zhu!dMq!iB-#;AL5a##R9KZn`{Gg(jpt<5;GC=)oRd{DHHyW=L%l>> zVu{h-vq=(daXE<&NavnSUo`GuY}`vFDkoB5cgX3JsSK8MOr|GT(gzuEuoBY$0p+o70E3NXC{e08B?_41%b+44C`K^Gh!)+4NVdc2`1#^ zl8TRS#o+PLbKk4V6T)EJ2j70dlD;TxkwoQEQkd^5Eg+;KMNBcyI|9pyO}amZiibN8 z5?BnfLHoYcJix`T;LixL)sK$d2Zq=p20c(5oB&UW{!hnl9fM9?KN3y3VyD1^4}r#E^%>Y>`s(6qDFr9?9ggk%uq24JOLVgTz;r-y>G&SCxO} zZ-%GN0j4SDHRQkZ7%$)cHs}`u5;mWDV~G6=n-ci=xPkQ zf-SB!)y*M)Q(;p}tABk#V@-9le|>4Fsiv^Hwy?&Z5__!xTlB$ffB%OT05tA!8{{)> zd-4Lx-#sYII)YKV=+IC8(pYF?kUUCi!fCPW%)H(4>i2)Jp<`?Kvh$HquYQoX{nsbm z4^7@1@!$QbyYsr@v|psH-E;KYAN_~7^e3VA&w9`&ii_TvpT>0o7a4cra^Wh%Rf=m7 zu9djP;KD!lwg~L~mN@fLKa_ovoO^eZWsq&c}|MtR_b02%=y|jJS&#OBQj(6QX z`^w>8NA~{v?w#|BHb0m;bFKZ#4#+->>m^*Lab3cd1f`>J(ML3ManTnZ{BtCb`@6W!!7DxhD-FPK@f{RpAa2E#R@$VJF-+IiK~jtw}}>aQ2R=e zFK=S84Ao!0XBK%nsK1Q*0iJi;@cMoDXBcb_FM{KZG1LLO~YNN*5VRLGIvyhynr^zMYf<`b4)x!w}> z&c>Z83vcRg07Q5TCIccgxh_Bza0DRodG=zYF4jgtIUr^tcN`${ba*`)Gr_M0B>t^} zo-ODe2th@`Q-IX&2Xun267(gkQyuiXfYe?Mh^FCN;RMj!5rBwwL56_IfFx%V^e1qV z+CLKTjDYl|I`Ll;5C4$fS8zE`t%0&hwl~eK|nmF zFZc*CCHSU*zZ39z0e>tYWirXF0wjM|0wQDu^jZU|au)(3=?mrpqU3%970FM5QM4zG z6Ywp}Gtnmnd{V$C1l%Uzg92_6aGijK0!|W;9vhIpR|NcL0lz0;w}4#&&KB^S0_F?o z67XKcmGrF@@Rx`m(f=;slL9^=;A#PD1-ui%Ci&Y1r1&6CNwWcwLP?VWseP<~HwgGs zK;%`0fR0^)7SOR%&;mMk2wFhLKM7hu$K!$)(D9g{1$1l^w1AF>1TCQBK|u@X*eYlN z9h(I$Ah(ZoQ$h9txqYI^KG9^CXaTu>qRBqBlUNTyRI3gwB81d-qTE&@8BnMD)myac)@fp#a*-3fGK0)1}+{lf(M zSOWb*0zETD(gQEW9fHa7;X*ULWq0sVS^poRGoUdeMM){A2wy)EZ-o4-0x+u*Bh zYH9Q1P53~W>Uw&i%{ifF{O=8`g?{p-7DwfP@Wls1b^bsA8tYr}hxWmG=JR(nLRsO} i17Y~f`GgnTLWSh>SJ!wGKH#mk*H=jY(}M*lh5ipq{f7wv literal 0 HcmV?d00001 diff --git a/src/macros.h b/src/macros.h new file mode 100644 index 0000000..9d98802 --- /dev/null +++ b/src/macros.h @@ -0,0 +1,22 @@ +#ifndef SEXP_MACROS +#define SEXP_MACROS + +#include +#include + +// define some macros +#define EXTRACT_NUMERIC_MATRIX(S_EXP, C_DATA, N_ROW, N_COL) \ + double *C_DATA = REAL(S_EXP); \ + const R_len_t N_ROW = nrows(S_EXP); \ + const R_len_t N_COL = ncols(S_EXP); + +#define EXTRACT_NUMERIC_VECTOR(S_EXP, C_DATA, LENGTH) \ + double *C_DATA = REAL(S_EXP); \ + const R_len_t LENGTH = length(S_EXP); + +#define ALLOC_VECTOR(type, size) (PROTECT(allocVector(type, size))) +#define ALLOC_REAL_VECTOR(size) (ALLOC_VECTOR(REALSXP, size)) +#define ALLOC_INTEGER_VECTOR(size) (ALLOC_VECTOR(INTSXP, size)) +#define ALLOC_LIST(size) (ALLOC_VECTOR(VECSXP, size)) + +#endif