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 0000000..baef9fb Binary files /dev/null and b/src/avl.o differ diff --git a/src/computeHypervolume.c b/src/computeHypervolume.c new file mode 100644 index 0000000..196f85c --- /dev/null +++ b/src/computeHypervolume.c @@ -0,0 +1,29 @@ +#include +#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 0000000..f51bdb1 Binary files /dev/null and b/src/computeHypervolume.o differ diff --git a/src/ecr.so b/src/ecr.so new file mode 100755 index 0000000..3845c97 Binary files /dev/null and b/src/ecr.so differ 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 0000000..188daf5 Binary files /dev/null and b/src/hv.o differ 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