PaGMO  1.1.5
fpl_cpp_original/hv.c
1 /*************************************************************************
2 
3  hypervolume computation
4 
5  ---------------------------------------------------------------------
6 
7  Copyright (c) 2010
8  Carlos M. Fonseca <cmfonsec@dei.uc.pt>
9  Manuel Lopez-Ibanez <manuel.lopez-ibanez@ulb.ac.be>
10  Luis Paquete <paquete@dei.uc.pt>
11  Andreia P. Guerreiro <andreia.guerreiro@ist.utl.pt>
12 
13  This program is free software (software libre); you can redistribute
14  it and/or modify it under the terms of the GNU General Public License
15  as published by the Free Software Foundation; either version 2 of the
16  License, or (at your option) any later version.
17 
18  This program is distributed in the hope that it will be useful, but
19  WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21  General Public License for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with this program; if not, you can obtain a copy of the GNU
25  General Public License at:
26  http://www.gnu.org/copyleft/gpl.html
27  or by writing to:
28  Free Software Foundation, Inc., 59 Temple Place,
29  Suite 330, Boston, MA 02111-1307 USA
30 
31  ----------------------------------------------------------------------
32 
33  Relevant literature:
34 
35  [1] C. M. Fonseca, L. Paquete, and M. Lopez-Ibanez. An
36  improved dimension-sweep algorithm for the hypervolume
37  indicator. In IEEE Congress on Evolutionary Computation,
38  pages 1157-1163, Vancouver, Canada, July 2006.
39 
40  [2] Nicola Beume, Carlos M. Fonseca, Manuel López-Ibáñez, Luís
41  Paquete, and J. Vahrenhold. On the complexity of computing the
42  hypervolume indicator. IEEE Transactions on Evolutionary
43  Computation, 13(5):1075-1082, 2009.
44 
45 *************************************************************************/
46 
47 #include "hv.h"
48 #include <stdlib.h>
49 #include <stdio.h>
50 #include <limits.h>
51 #include <float.h>
52 #include <assert.h>
53 
54 
55 static int compare_tree_asc(const void *p1, const void *p2);
56 
57 /*-----------------------------------------------------------------------------
58 
59  The following is a reduced version of the AVL-tree library used here
60  according to the terms of the GPL. See the copyright notice below.
61 
62 */
63 #define AVL_DEPTH
64 
65 /*****************************************************************************
66 
67  avl.h - Source code for the AVL-tree library.
68 
69  Copyright (C) 1998 Michael H. Buselli <cosine@cosine.org>
70  Copyright (C) 2000-2002 Wessel Dankers <wsl@nl.linux.org>
71 
72  This library is free software; you can redistribute it and/or
73  modify it under the terms of the GNU Lesser General Public
74  License as published by the Free Software Foundation; either
75  version 2.1 of the License, or (at your option) any later version.
76 
77  This library is distributed in the hope that it will be useful,
78  but WITHOUT ANY WARRANTY; without even the implied warranty of
79  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
80  Lesser General Public License for more details.
81 
82  You should have received a copy of the GNU Lesser General Public
83  License along with this library; if not, write to the Free Software
84  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
85 
86  Augmented AVL-tree. Original by Michael H. Buselli <cosine@cosine.org>.
87 
88  Modified by Wessel Dankers <wsl@nl.linux.org> to add a bunch of bloat to
89  the sourcecode, change the interface and squash a few bugs.
90  Mail him if you find new bugs.
91 
92 *****************************************************************************/
93 
94 /* User supplied function to compare two items like strcmp() does.
95  * For example: cmp(a,b) will return:
96  * -1 if a < b
97  * 0 if a = b
98  * 1 if a > b
99  */
100 typedef int (*avl_compare_t)(const void *, const void *);
101 
102 /* User supplied function to delete an item when a node is free()d.
103  * If NULL, the item is not free()d.
104  */
105 typedef void (*avl_freeitem_t)(void *);
106 
107 typedef struct avl_node_t {
108  struct avl_node_t *next;
109  struct avl_node_t *prev;
110  struct avl_node_t *parent;
111  struct avl_node_t *left;
112  struct avl_node_t *right;
113  void *item;
114  double domr;
115 #ifdef AVL_DEPTH
116  unsigned char depth;
117 #endif
118 } avl_node_t;
119 
120 typedef struct avl_tree_t {
121  avl_node_t *head;
122  avl_node_t *tail;
123  avl_node_t *top;
124  avl_compare_t cmp;
125  avl_freeitem_t freeitem;
126 } avl_tree_t;
127 
128 
129 /*****************************************************************************
130 
131  avl.c - Source code for the AVL-tree library.
132 
133 *****************************************************************************/
134 
135 static void avl_rebalance(avl_tree_t *, avl_node_t *);
136 
137 #ifdef AVL_DEPTH
138 #define NODE_DEPTH(n) ((n) ? (n)->depth : 0)
139 #define L_DEPTH(n) (NODE_DEPTH((n)->left))
140 #define R_DEPTH(n) (NODE_DEPTH((n)->right))
141 #define CALC_DEPTH(n) ((L_DEPTH(n)>R_DEPTH(n)?L_DEPTH(n):R_DEPTH(n)) + 1)
142 #endif
143 
144 static int avl_check_balance(avl_node_t *avlnode) {
145 #ifdef AVL_DEPTH
146  int d;
147  d = R_DEPTH(avlnode) - L_DEPTH(avlnode);
148  return d<-1?-1:d>1?1:0;
149 #endif
150 }
151 
152 static int
153 avl_search_closest(const avl_tree_t *avltree, const void *item, avl_node_t **avlnode) {
154  avl_node_t *node;
155  int c;
156 
157  if(!avlnode)
158  avlnode = &node;
159 
160  node = avltree->top;
161 
162  if(!node)
163  return *avlnode = NULL, 0;
164 
165 
166  for(;;) {
167  c = compare_tree_asc(item, node->item);
168 
169  if(c < 0) {
170  if(node->left)
171  node = node->left;
172  else
173  return *avlnode = node, -1;
174  } else if(c > 0) {
175  if(node->right)
176  node = node->right;
177  else
178  return *avlnode = node, 1;
179  } else {
180  return *avlnode = node, 0;
181  }
182  }
183 }
184 
185 static avl_tree_t *
186 avl_init_tree(avl_tree_t *rc, avl_compare_t cmp, avl_freeitem_t freeitem) {
187  if(rc) {
188  rc->head = NULL;
189  rc->tail = NULL;
190  rc->top = NULL;
191  rc->cmp = cmp;
192  rc->freeitem = freeitem;
193  }
194  return rc;
195 }
196 
197 static avl_tree_t *
198 avl_alloc_tree(avl_compare_t cmp, avl_freeitem_t freeitem) {
199  return avl_init_tree(malloc(sizeof(avl_tree_t)), cmp, freeitem);
200 }
201 
202 static void
203 avl_clear_tree(avl_tree_t *avltree) {
204  avltree->top = avltree->head = avltree->tail = NULL;
205 }
206 
207 static void
208 avl_clear_node(avl_node_t *newnode) {
209  newnode->left = newnode->right = NULL;
210  #ifdef AVL_COUNT
211  newnode->count = 1;
212  #endif
213  #ifdef AVL_DEPTH
214  newnode->depth = 1;
215  #endif
216 }
217 
218 static avl_node_t *
219 avl_insert_top(avl_tree_t *avltree, avl_node_t *newnode) {
220  avl_clear_node(newnode);
221  newnode->prev = newnode->next = newnode->parent = NULL;
222  avltree->head = avltree->tail = avltree->top = newnode;
223  return newnode;
224 }
225 
226 static avl_node_t *
227 avl_insert_before(avl_tree_t *avltree, avl_node_t *node, avl_node_t *newnode) {
228 /* if(!node)
229  return avltree->tail
230  ? avl_insert_after(avltree, avltree->tail, newnode)
231  : avl_insert_top(avltree, newnode);
232 
233  if(node->left)
234  return avl_insert_after(avltree, node->prev, newnode);
235 */
236  assert (node);
237  assert (!node->left);
238 
239  avl_clear_node(newnode);
240 
241  newnode->next = node;
242  newnode->parent = node;
243 
244  newnode->prev = node->prev;
245  if(node->prev)
246  node->prev->next = newnode;
247  else
248  avltree->head = newnode;
249  node->prev = newnode;
250 
251  node->left = newnode;
252  avl_rebalance(avltree, node);
253  return newnode;
254 }
255 
256 static avl_node_t *
257 avl_insert_after(avl_tree_t *avltree, avl_node_t *node, avl_node_t *newnode) {
258 /* if(!node)
259  return avltree->head
260  ? avl_insert_before(avltree, avltree->head, newnode)
261  : avl_insert_top(avltree, newnode);
262 
263  if(node->right)
264  return avl_insert_before(avltree, node->next, newnode);
265 */
266  assert (node);
267  assert (!node->right);
268 
269  avl_clear_node(newnode);
270 
271  newnode->prev = node;
272  newnode->parent = node;
273 
274  newnode->next = node->next;
275  if(node->next)
276  node->next->prev = newnode;
277  else
278  avltree->tail = newnode;
279  node->next = newnode;
280 
281  node->right = newnode;
282  avl_rebalance(avltree, node);
283  return newnode;
284 }
285 
286 /*
287  * avl_unlink_node:
288  * Removes the given node. Does not delete the item at that node.
289  * The item of the node may be freed before calling avl_unlink_node.
290  * (In other words, it is not referenced by this function.)
291  */
292 static void
293 avl_unlink_node(avl_tree_t *avltree, avl_node_t *avlnode) {
294  avl_node_t *parent;
295  avl_node_t **superparent;
296  avl_node_t *subst, *left, *right;
297  avl_node_t *balnode;
298 
299  if(avlnode->prev)
300  avlnode->prev->next = avlnode->next;
301  else
302  avltree->head = avlnode->next;
303 
304  if(avlnode->next)
305  avlnode->next->prev = avlnode->prev;
306  else
307  avltree->tail = avlnode->prev;
308 
309  parent = avlnode->parent;
310 
311  superparent = parent
312  ? avlnode == parent->left ? &parent->left : &parent->right
313  : &avltree->top;
314 
315  left = avlnode->left;
316  right = avlnode->right;
317  if(!left) {
318  *superparent = right;
319  if(right)
320  right->parent = parent;
321  balnode = parent;
322  } else if(!right) {
323  *superparent = left;
324  left->parent = parent;
325  balnode = parent;
326  } else {
327  subst = avlnode->prev;
328  if(subst == left) {
329  balnode = subst;
330  } else {
331  balnode = subst->parent;
332  balnode->right = subst->left;
333  if(balnode->right)
334  balnode->right->parent = balnode;
335  subst->left = left;
336  left->parent = subst;
337  }
338  subst->right = right;
339  subst->parent = parent;
340  right->parent = subst;
341  *superparent = subst;
342  }
343 
344  avl_rebalance(avltree, balnode);
345 }
346 
347 /*
348  * avl_rebalance:
349  * Rebalances the tree if one side becomes too heavy. This function
350  * assumes that both subtrees are AVL-trees with consistant data. The
351  * function has the additional side effect of recalculating the count of
352  * the tree at this node. It should be noted that at the return of this
353  * function, if a rebalance takes place, the top of this subtree is no
354  * longer going to be the same node.
355  */
356 static void
357 avl_rebalance(avl_tree_t *avltree, avl_node_t *avlnode) {
358  avl_node_t *child;
359  avl_node_t *gchild;
360  avl_node_t *parent;
361  avl_node_t **superparent;
362 
363  parent = avlnode;
364 
365  while(avlnode) {
366  parent = avlnode->parent;
367 
368  superparent = parent
369  ? avlnode == parent->left ? &parent->left : &parent->right
370  : &avltree->top;
371 
372  switch(avl_check_balance(avlnode)) {
373  case -1:
374  child = avlnode->left;
375  #ifdef AVL_DEPTH
376  if(L_DEPTH(child) >= R_DEPTH(child)) {
377  #else
378  #ifdef AVL_COUNT
379  if(L_COUNT(child) >= R_COUNT(child)) {
380  #else
381  #error No balancing possible.
382  #endif
383  #endif
384  avlnode->left = child->right;
385  if(avlnode->left)
386  avlnode->left->parent = avlnode;
387  child->right = avlnode;
388  avlnode->parent = child;
389  *superparent = child;
390  child->parent = parent;
391  #ifdef AVL_COUNT
392  avlnode->count = CALC_COUNT(avlnode);
393  child->count = CALC_COUNT(child);
394  #endif
395  #ifdef AVL_DEPTH
396  avlnode->depth = CALC_DEPTH(avlnode);
397  child->depth = CALC_DEPTH(child);
398  #endif
399  } else {
400  gchild = child->right;
401  avlnode->left = gchild->right;
402  if(avlnode->left)
403  avlnode->left->parent = avlnode;
404  child->right = gchild->left;
405  if(child->right)
406  child->right->parent = child;
407  gchild->right = avlnode;
408  if(gchild->right)
409  gchild->right->parent = gchild;
410  gchild->left = child;
411  if(gchild->left)
412  gchild->left->parent = gchild;
413  *superparent = gchild;
414  gchild->parent = parent;
415  #ifdef AVL_COUNT
416  avlnode->count = CALC_COUNT(avlnode);
417  child->count = CALC_COUNT(child);
418  gchild->count = CALC_COUNT(gchild);
419  #endif
420  #ifdef AVL_DEPTH
421  avlnode->depth = CALC_DEPTH(avlnode);
422  child->depth = CALC_DEPTH(child);
423  gchild->depth = CALC_DEPTH(gchild);
424  #endif
425  }
426  break;
427  case 1:
428  child = avlnode->right;
429  #ifdef AVL_DEPTH
430  if(R_DEPTH(child) >= L_DEPTH(child)) {
431  #else
432  #ifdef AVL_COUNT
433  if(R_COUNT(child) >= L_COUNT(child)) {
434  #else
435  #error No balancing possible.
436  #endif
437  #endif
438  avlnode->right = child->left;
439  if(avlnode->right)
440  avlnode->right->parent = avlnode;
441  child->left = avlnode;
442  avlnode->parent = child;
443  *superparent = child;
444  child->parent = parent;
445  #ifdef AVL_COUNT
446  avlnode->count = CALC_COUNT(avlnode);
447  child->count = CALC_COUNT(child);
448  #endif
449  #ifdef AVL_DEPTH
450  avlnode->depth = CALC_DEPTH(avlnode);
451  child->depth = CALC_DEPTH(child);
452  #endif
453  } else {
454  gchild = child->left;
455  avlnode->right = gchild->left;
456  if(avlnode->right)
457  avlnode->right->parent = avlnode;
458  child->left = gchild->right;
459  if(child->left)
460  child->left->parent = child;
461  gchild->left = avlnode;
462  if(gchild->left)
463  gchild->left->parent = gchild;
464  gchild->right = child;
465  if(gchild->right)
466  gchild->right->parent = gchild;
467  *superparent = gchild;
468  gchild->parent = parent;
469  #ifdef AVL_COUNT
470  avlnode->count = CALC_COUNT(avlnode);
471  child->count = CALC_COUNT(child);
472  gchild->count = CALC_COUNT(gchild);
473  #endif
474  #ifdef AVL_DEPTH
475  avlnode->depth = CALC_DEPTH(avlnode);
476  child->depth = CALC_DEPTH(child);
477  gchild->depth = CALC_DEPTH(gchild);
478  #endif
479  }
480  break;
481  default:
482  #ifdef AVL_COUNT
483  avlnode->count = CALC_COUNT(avlnode);
484  #endif
485  #ifdef AVL_DEPTH
486  avlnode->depth = CALC_DEPTH(avlnode);
487  #endif
488  }
489  avlnode = parent;
490  }
491 }
492 
493 /*------------------------------------------------------------------------------
494  end of functions from AVL-tree library.
495 *******************************************************************************/
496 #define VARIANT 4
497 #if !defined(VARIANT) || VARIANT < 1 || VARIANT > 4
498 #error VARIANT must be either 1, 2, 3 or 4, e.g., 'make VARIANT=4'
499 #endif
500 
501 #if __GNUC__ >= 3
502 # define __hv_unused __attribute__ ((unused))
503 #else
504 # define __hv_unused /* no 'unused' attribute available */
505 #endif
506 
507 #if VARIANT < 3
508 # define __variant3_only __hv_unused
509 #else
510 # define __variant3_only
511 #endif
512 
513 #if VARIANT < 2
514 # define __variant2_only __hv_unused
515 #else
516 # define __variant2_only
517 #endif
518 
519 typedef struct dlnode {
520  double *x; /* The data vector */
521  struct dlnode **next; /* Next-node vector */
522  struct dlnode **prev; /* Previous-node vector */
523  struct avl_node_t * tnode;
524  int ignore;
525  int ignore_best; //used in define_order
526 #if VARIANT >= 2
527  double *area; /* Area */
528 #endif
529 #if VARIANT >= 3
530  double *vol; /* Volume */
531 #endif
532 } dlnode_t;
533 
534 static avl_tree_t *tree;
535 #if VARIANT < 4
536 int stop_dimension = 1; /* default: stop on dimension 2 */
537 #else
538 int stop_dimension = 2; /* default: stop on dimension 3 */
539 #endif
540 
541 static int compare_node(const void *p1, const void* p2)
542 {
543  const double x1 = *((*(const dlnode_t **)p1)->x);
544  const double x2 = *((*(const dlnode_t **)p2)->x);
545 
546  return (x1 < x2) ? -1 : (x1 > x2) ? 1 : 0;
547 }
548 
549 static int compare_tree_asc(const void *p1, const void *p2)
550 {
551  const double *x1 = (const double *)p1;
552  const double *x2 = (const double *)p2;
553 
554  return (x1[1] > x2[1]) ? -1 : (x1[1] < x2[1]) ? 1
555  : (x1[0] >= x2[0]) ? -1 : 1;
556 }
557 
558 /*
559  * Setup circular double-linked list in each dimension
560  */
561 
562 static dlnode_t *
563 setup_cdllist(double *data, int d, int n)
564 {
565  dlnode_t *head;
566  dlnode_t **scratch;
567  int i, j;
568 
569  head = malloc ((n+1) * sizeof(dlnode_t));
570 
571  head->x = data;
572  head->ignore = 0; /* should never get used */
573  head->next = malloc( d * (n+1) * sizeof(dlnode_t*));
574  head->prev = malloc( d * (n+1) * sizeof(dlnode_t*));
575  head->tnode = malloc ((n+1) * sizeof(avl_node_t));
576 
577 #if VARIANT >= 2
578  head->area = malloc(d * (n+1) * sizeof(double));
579 #endif
580 #if VARIANT >= 3
581  head->vol = malloc(d * (n+1) * sizeof(double));
582 #endif
583 
584  for (i = 1; i <= n; i++) {
585  head[i].x = head[i-1].x + d;/* this will be fixed a few lines below... */
586  head[i].ignore = 0;
587  head[i].next = head[i-1].next + d;
588  head[i].prev = head[i-1].prev + d;
589  head[i].tnode = head[i-1].tnode + 1;
590 #if VARIANT >= 2
591  head[i].area = head[i-1].area + d;
592 #endif
593 #if VARIANT >= 3
594  head[i].vol = head[i-1].vol + d;
595 #endif
596  }
597  head->x = NULL; /* head contains no data */
598 
599  scratch = malloc(n * sizeof(dlnode_t*));
600  for (i = 0; i < n; i++)
601  scratch[i] = head + i + 1;
602 
603  for (j = d-1; j >= 0; j--) {
604  for (i = 0; i < n; i++)
605  scratch[i]->x--;
606  qsort(scratch, n, sizeof(dlnode_t*), compare_node);
607  head->next[j] = scratch[0];
608  scratch[0]->prev[j] = head;
609  for (i = 1; i < n; i++) {
610  scratch[i-1]->next[j] = scratch[i];
611  scratch[i]->prev[j] = scratch[i-1];
612  }
613  scratch[n-1]->next[j] = head;
614  head->prev[j] = scratch[n-1];
615  }
616 
617  free(scratch);
618 
619  for (i = 1; i <= n; i++) {
620  (head[i].tnode)->item = head[i].x;
621  }
622 
623  return head;
624 }
625 
626 static void free_cdllist(dlnode_t * head)
627 {
628  free(head->tnode); /* Frees _all_ nodes. */
629  free(head->next);
630  free(head->prev);
631 #if VARIANT >= 2
632  free(head->area);
633 #endif
634 #if VARIANT >= 3
635  free(head->vol);
636 #endif
637  free(head);
638 }
639 
640 static void delete (dlnode_t *nodep, int dim, double * bound __variant3_only)
641 {
642  int i;
643 
644  for (i = stop_dimension; i < dim; i++) {
645  nodep->prev[i]->next[i] = nodep->next[i];
646  nodep->next[i]->prev[i] = nodep->prev[i];
647 #if VARIANT >= 3
648  if (bound[i] > nodep->x[i])
649  bound[i] = nodep->x[i];
650 #endif
651  }
652 }
653 
654 #if VARIANT >= 2
655 static void delete_dom (dlnode_t *nodep, int dim)
656 {
657  int i;
658 
659  for (i = stop_dimension; i < dim; i++) {
660  nodep->prev[i]->next[i] = nodep->next[i];
661  nodep->next[i]->prev[i] = nodep->prev[i];
662  }
663 }
664 #endif
665 
666 static void reinsert (dlnode_t *nodep, int dim, double * bound __variant3_only)
667 {
668  int i;
669 
670  for (i = stop_dimension; i < dim; i++) {
671  nodep->prev[i]->next[i] = nodep;
672  nodep->next[i]->prev[i] = nodep;
673 #if VARIANT >= 3
674  if (bound[i] > nodep->x[i])
675  bound[i] = nodep->x[i];
676 #endif
677  }
678 }
679 
680 #if VARIANT >= 2
681 static void reinsert_dom (dlnode_t *nodep, int dim)
682 {
683  int i;
684  for (i = stop_dimension; i < dim; i++) {
685  dlnode_t *p = nodep->prev[i];
686  p->next[i] = nodep;
687  nodep->next[i]->prev[i] = nodep;
688  nodep->area[i] = p->area[i];
689 
690  #if VARIANT >= 3
691  nodep->vol[i] = p->vol[i] + p->area[i] * (nodep->x[i] - p->x[i]);
692  #endif
693  }
694 }
695 #endif
696 
697 static double
698 hv_recursive(dlnode_t *list, int dim, int c, const double * ref,
699  double * bound)
700 {
701  /* ------------------------------------------------------
702  General case for dimensions higher than stop_dimension
703  ------------------------------------------------------ */
704  if ( dim > stop_dimension ) {
705  dlnode_t *p0 = list;
706  dlnode_t *p1 = list->prev[dim];
707  double hyperv = 0;
708 #if VARIANT == 1
709  double hypera;
710 #endif
711 #if VARIANT >= 2
712  dlnode_t *pp;
713  for (pp = p1; pp->x; pp = pp->prev[dim]) {
714  if (pp->ignore < dim)
715  pp->ignore = 0;
716  }
717 #endif
718  while (c > 1
719 #if VARIANT >= 3
720  /* We delete all points x[dim] > bound[dim]. In case of
721  repeated coordinates, we also delete all points
722  x[dim] == bound[dim] except one. */
723  && (p1->x[dim] > bound[dim]
724  || p1->prev[dim]->x[dim] >= bound[dim])
725 #endif
726  ) {
727  p0 = p1;
728 #if VARIANT >=2
729  if (p0->ignore >= dim)
730  delete_dom(p0, dim);
731  else
732  delete(p0, dim, bound);
733 #else
734  delete(p0, dim, bound);
735 #endif
736  p1 = p0->prev[dim];
737  c--;
738  }
739 
740 #if VARIANT == 1
741  hypera = hv_recursive(list, dim-1, c, ref, bound);
742 
743 #elif VARIANT == 2
744  int i;
745  p1->area[0] = 1;
746  for (i = 1; i <= dim; i++)
747  p1->area[i] = p1->area[i-1] * (ref[i-1] - p1->x[i-1]);
748 
749 #elif VARIANT >= 3
750  if (c > 1) {
751  hyperv = p1->prev[dim]->vol[dim] + p1->prev[dim]->area[dim]
752  * (p1->x[dim] - p1->prev[dim]->x[dim]);
753 
754  if (p1->ignore >= dim)
755  p1->area[dim] = p1->prev[dim]->area[dim];
756  else {
757  p1->area[dim] = hv_recursive(list, dim - 1, c, ref, bound);
758  /* At this point, p1 is the point with the highest value in
759  dimension dim in the list, so if it is dominated in
760  dimension dim-1, so it is also dominated in dimension
761  dim. */
762  if (p1->ignore == (dim - 1))
763  p1->ignore = dim;
764  }
765  } else {
766  int i;
767  p1->area[0] = 1;
768  for (i = 1; i <= dim; i++)
769  p1->area[i] = p1->area[i-1] * (ref[i-1] - p1->x[i-1]);
770  }
771  p1->vol[dim] = hyperv;
772 #endif
773 
774  while (p0->x != NULL) {
775 
776 #if VARIANT == 1
777  hyperv += hypera * (p0->x[dim] - p1->x[dim]);
778 #else
779  hyperv += p1->area[dim] * (p0->x[dim] - p1->x[dim]);
780 #endif
781  c++;
782 #if VARIANT >= 2
783  if (p0->ignore >= dim) {
784  reinsert_dom (p0, dim);
785  p0->area[dim] = p1->area[dim];
786  } else {
787 #endif
788  reinsert (p0, dim, bound);
789 #if VARIANT >= 2
790  p0->area[dim] = hv_recursive (list, dim-1, c, ref, bound);
791  if (p0->ignore == (dim - 1))
792  p0->ignore = dim;
793  }
794 #elif VARIANT == 1
795  hypera = hv_recursive (list, dim-1, c, ref, NULL);
796 #endif
797  p1 = p0;
798  p0 = p0->next[dim];
799 #if VARIANT >= 3
800  p1->vol[dim] = hyperv;
801 #endif
802  }
803 #if VARIANT >= 3
804  bound[dim] = p1->x[dim];
805 #endif
806 
807 #if VARIANT == 1
808  hyperv += hypera * (ref[dim] - p1->x[dim]);
809 #else
810  hyperv += p1->area[dim] * (ref[dim] - p1->x[dim]);
811 #endif
812  return hyperv;
813  }
814 
815 
816  /* ---------------------------
817  special case of dimension 3
818  --------------------------- */
819  else if (dim == 2) {
820  double hyperv;
821  double hypera;
822  double height;
823 
824 #if VARIANT >= 3
825  dlnode_t *pp = list->prev[2];
826  avl_node_t *tnode;
827 
828  /* All the point that have value of x[2] lower than bound[2] are points
829  that were previously processed, so there's no need to process them
830  again. In this case, every point was processed before, so the
831  volume is known. */
832  if (pp->x[2] < bound[2])
833  return pp->vol[2] + pp->area[2] * (ref[2] - pp->x[2]);
834 
835  pp = list->next[2];
836 
837  /* In this case, every point has to be processed. */
838  if (pp->x[2] >= bound[2]) {
839  pp->tnode->domr = ref[2];
840  pp->area[2] = (ref[0] - pp->x[0]) * (ref[1] - pp->x[1]);
841  pp->vol[2] = 0;
842  pp->ignore = 0;
843  } else {
844  /* Otherwise, we look for the first point that has to be in the
845  tree, by searching for the first point that isn't dominated or
846  that is dominated by a point with value of x[2] higher or equal
847  than bound[2] (domr keeps the value of the x[2] of the point
848  that dominates pp, or ref[2] if it isn't dominated). */
849  while (pp->tnode->domr < bound[2]) {
850  pp = pp->next[2];
851  }
852  }
853 
854  pp->ignore = 0;
855  avl_insert_top(tree,pp->tnode);
856  pp->tnode->domr = ref[2];
857 
858  /* Connect all points that aren't dominated or that are dominated and
859  the point that dominates it has value x[2] (pp->tnode->domr) equal
860  or higher than bound[2]. */
861  for (pp = pp->next[2]; pp->x[2] < bound[2]; pp = pp->next[2]) {
862  if (pp->tnode->domr >= bound[2]) {
863  avl_node_t *tnodeaux = pp->tnode;
864  tnodeaux->domr = ref[2];
865  if (avl_search_closest(tree, pp->x, &tnode) <= 0)
866  avl_insert_before(tree, tnode, tnodeaux);
867  else
868  avl_insert_after(tree, tnode, tnodeaux);
869  }
870  }
871  pp = pp->prev[2];
872  hyperv = pp->vol[2];
873  hypera = pp->area[2];
874 
875  height = (pp->next[2]->x)
876  ? pp->next[2]->x[2] - pp->x[2]
877  : ref[2] - pp->x[2];
878 
879  bound[2] = list->prev[2]->x[2];
880 #else
881 /* VARIANT <= 2 */
882  dlnode_t *pp = list->next[2];
883  hyperv = 0;
884 
885  hypera = (ref[0] - pp->x[0])*(ref[1] - pp->x[1]);
886  height = (c == 1)
887  ? ref[2] - pp->x[2]
888  : pp->next[2]->x[2] - pp->x[2];
889 
890  avl_insert_top(tree,pp->tnode);
891 #endif
892  hyperv += hypera * height;
893  for (pp = pp->next[2]; pp->x != NULL; pp = pp->next[2]) {
894  const double * prv_ip, * nxt_ip;
895  avl_node_t *tnode;
896  int cmp;
897 
898 #if VARIANT >= 3
899  pp->vol[2] = hyperv;
900 #endif
901  height = (pp == list->prev[2])
902  ? ref[2] - pp->x[2]
903  : pp->next[2]->x[2] - pp->x[2];
904 #if VARIANT >= 2
905  if (pp->ignore >= 2) {
906  hyperv += hypera * height;
907 #if VARIANT >= 3
908  pp->area[2] = hypera;
909 #endif
910  continue;
911  }
912 #endif
913  cmp = avl_search_closest(tree, pp->x, &tnode);
914  if (cmp <= 0) {
915  nxt_ip = (double *)(tnode->item);
916  } else {
917  nxt_ip = (tnode->next != NULL)
918  ? (double *)(tnode->next->item)
919  : ref;
920  }
921  if (nxt_ip[0] <= pp->x[0]) {
922  pp->ignore = 2;
923 #if VARIANT >= 3
924  pp->tnode->domr = pp->x[2];
925  pp->area[2] = hypera;
926 #endif
927  if (height > 0)
928  hyperv += hypera * height;
929  continue;
930  }
931  if (cmp <= 0) {
932  avl_insert_before(tree, tnode, pp->tnode);
933  tnode = pp->tnode->prev;
934  } else {
935  avl_insert_after(tree, tnode, pp->tnode);
936  }
937 #if VARIANT >= 3
938  pp->tnode->domr = ref[2];
939 #endif
940  if (tnode != NULL) {
941  prv_ip = (double *)(tnode->item);
942  if (prv_ip[0] >= pp->x[0]) {
943  const double * cur_ip;
944 
945  tnode = pp->tnode->prev;
946  /* cur_ip = point dominated by pp with highest
947  [0]-coordinate. */
948  cur_ip = (double *)(tnode->item);
949  while (tnode->prev) {
950  prv_ip = (double *)(tnode->prev->item);
951  hypera -= (prv_ip[1] - cur_ip[1]) * (nxt_ip[0] - cur_ip[0]);
952  if (prv_ip[0] < pp->x[0])
953  break; /* prv is not dominated by pp */
954  cur_ip = prv_ip;
955  avl_unlink_node(tree,tnode);
956 #if VARIANT >= 3
957  /* saves the value of x[2] of the point that
958  dominates tnode. */
959  tnode->domr = pp->x[2];
960 #endif
961  tnode = tnode->prev;
962  }
963 
964  avl_unlink_node(tree, tnode);
965 #if VARIANT >= 3
966  tnode->domr = pp->x[2];
967 #endif
968  if (!tnode->prev) {
969  hypera -= (ref[1] - cur_ip[1]) * (nxt_ip[0] - cur_ip[0]);
970  prv_ip = ref;
971  }
972  }
973  } else
974  prv_ip = ref;
975 
976  hypera += (prv_ip[1] - pp->x[1]) * (nxt_ip[0] - pp->x[0]);
977 
978  if (height > 0)
979  hyperv += hypera * height;
980 #if VARIANT >= 3
981  pp->area[2] = hypera;
982 #endif
983  }
984  avl_clear_tree(tree);
985  return hyperv;
986  }
987 
988  /* special case of dimension 2 */
989  else if (dim == 1) {
990  const dlnode_t *p1 = list->next[1];
991  double hypera = p1->x[0];
992  double hyperv = 0;
993  dlnode_t *p0;
994 
995  while ((p0 = p1->next[1])->x) {
996  hyperv += (ref[0] - hypera) * (p0->x[1] - p1->x[1]);
997  if (p0->x[0] < hypera)
998  hypera = p0->x[0];
999  else if (p0->ignore == 0)
1000  p0->ignore = 1;
1001  p1 = p0;
1002  }
1003  hyperv += (ref[0] - hypera) * (ref[1] - p1->x[1]);
1004  return hyperv;
1005  }
1006 
1007  /* special case of dimension 1 */
1008  else if (dim == 0) {
1009  list->next[0]->ignore = -1;
1010  return (ref[0] - list->next[0]->x[0]);
1011  }
1012  else {
1013  fprintf(stderr, "%s:%d: unreachable condition! \n"
1014  "This is a bug, please report it to "
1015  "manuel.lopez-ibanez@ulb.ac.be\n", __FILE__, __LINE__);
1016  exit(EXIT_FAILURE);
1017  }
1018 }
1019 
1020 /*
1021  Removes the point from the circular double-linked list, but it
1022  doesn't remove the data.
1023 */
1024 static void
1025 filter_delete_node(dlnode_t *node, int d)
1026 {
1027  int i;
1028 
1029  for (i = 0; i < d; i++) {
1030  node->next[i]->prev[i] = node->prev[i];
1031  node->prev[i]->next[i] = node->next[i];
1032  }
1033 }
1034 
1035 /*
1036  Filters those points that do not strictly dominate the reference
1037  point. This is needed to assure that the points left are only those
1038  that are needed to calculate the hypervolume.
1039 */
1040 static int
1041 filter(dlnode_t *list, int d, int n, const double *ref)
1042 {
1043  int i, j;
1044 
1045  /* fprintf (stderr, "%d points initially\n", n); */
1046  for (i = 0; i < d; i++) {
1047  dlnode_t *aux = list->prev[i];
1048  int np = n;
1049  for (j = 0; j < np; j++) {
1050  if (aux->x[i] < ref[i])
1051  break;
1052  filter_delete_node (aux, d);
1053  aux = aux->prev[i];
1054  n--;
1055  }
1056  }
1057  /* fprintf (stderr, "%d points remain\n", n); */
1058  return n;
1059 }
1060 
1061 
1062 #ifdef EXPERIMENTAL
1063 /*
1064  Verifies up to which dimension k, domr dominates p and returns k
1065  (it is assumed that domr doesn't dominate p in dimensions higher than dim).
1066 */
1067 static int
1068 test_domr(dlnode_t *p, dlnode_t *domr, int dim, int *order)
1069 {
1070  int i;
1071  for(i = 1; i <= dim; i++){
1072  if (p->x[order[i]] < domr->x[order[i]])
1073  return i - 1;
1074  }
1075  return dim;
1076 }
1077 
1078 /*
1079  Verifies up to which dimension k the point pp is dominated and
1080  returns k. This functions is called only to verify points that
1081  aren't dominated for more than dim dimensions, so k will always be
1082  lower or equal to dim.
1083 */
1084 static int
1085 test_dom(dlnode_t *list, dlnode_t *pp, int dim, int *order)
1086 {
1087  dlnode_t *p0;
1088  int r, r_b = 0;
1089  int i = order[0];
1090 
1091  p0 = list->next[i];
1092 
1093  /* In every iteration, it is verified if p0 dominates pp and
1094  up to which dimension. The goal is to find the point that
1095  dominates pp in more dimension, starting in dimension 0.
1096 
1097  Points are processed in ascending order of the first
1098  dimension. This means that if a point p0 is dominated in
1099  the first k dimensions, where k >=dim, then the point that
1100  dominates it (in the first k dimensions) was already
1101  processed, so p0 won't dominate pp in more dimensions that
1102  the point that dominates p0 (because pp can be dominated,
1103  at most, up to dim dimensions, and so if p0 dominates pp in
1104  the first y dimensions (y < dim), the point that dominates
1105  p0 also dominates pp in the first y dimensions or more, and
1106  this informations is already stored in r_b), so p0 is
1107  skipped. */
1108  while (p0 != pp) {
1109  if (p0->ignore < dim) {
1110  r = test_domr (pp, p0, dim, order);
1111  /* if pp is dominated in the first dim + 1 dimensions,
1112  it is not necessary to verify other points that
1113  might dominate pp, because pp won't be dominated in
1114  more that dim+1 dimensions. */
1115  if (r == dim) return r;
1116  else if (r > r_b) r_b = r;
1117  }
1118  p0 = p0->next[i];
1119  }
1120  return r_b;
1121 }
1122 
1123 /*
1124  Determines the number of dominated points from dimension 0 to k,
1125  where k <= dim.
1126 */
1127 static void determine_ndom(dlnode_t *list, int dim, int *order, int *count)
1128 {
1129  dlnode_t *p1;
1130  int i, dom;
1131  int ord = order[0];
1132 
1133  for (i = 0; i <= dim; i++)
1134  count[i] = 0;
1135 
1136  p1 = list->next[ord];
1137  p1->ignore = 0;
1138 
1139  p1 = list->next[ord];
1140 
1141  while (p1 != list) {
1142  if (p1->ignore <= dim) {
1143  dom = test_dom(list, p1, dim, order);
1144  count[dom]++;
1145  p1->ignore = dom;
1146  }
1147  p1 = p1->next[ord];
1148  }
1149 }
1150 
1151 static void delete_dominated(dlnode_t *nodep, int dim)
1152 {
1153  int i;
1154  for (i = 0; i <= dim; i++) {
1155  nodep->prev[i]->next[i] = nodep->next[i];
1156  nodep->next[i]->prev[i] = nodep->prev[i];
1157  }
1158 }
1159 
1160 
1161 /*
1162  Determines the number of dominated points from dimension 0 to k,
1163  where k <= dim, for the original order of objectives. Also defines
1164  that this order is the best order so far, so every point has the
1165  information up to which dimension it is dominated (ignore) and it is
1166  considered the highest number of dimensions in which it is dominated
1167  (so ignore_best is also updated).
1168 
1169  If there is any point dominated in every dimension, seen that it
1170  doesn't contribute to the hypervolume, it is removed as soon as
1171  possible, this way there's no waste of time with these points.
1172  Returns the number of total points. */
1173 static int
1174 determine_ndomf(dlnode_t *list, int dim, int c, int *order, int *count)
1175 {
1176  dlnode_t *p1;
1177  int i, dom;
1178  int ord = order[0];
1179 
1180  for(i = 0; i <= dim; i++)
1181  count[i] = 0;
1182 
1183  p1 = list->next[ord];
1184  p1->ignore = p1->ignore_best = 0;
1185 
1186  p1 = list->next[ord];
1187 
1188  /* Determines up to which dimension each point is dominated and
1189  uses this information to count the number of dominated points
1190  from dimension 0 to k, where k <= dim.
1191 
1192  Points that are dominated in more than the first 'dim'
1193  dimensions will continue to be dominated in those dimensions,
1194  and so they're skipped, it's not necessary to find out again up
1195  to which dimension they're dominated. */
1196  while (p1 != list){
1197  if (p1->ignore <= dim) {
1198  dom = test_dom(list, p1, dim, order);
1199  count[dom]++;
1200  p1->ignore = p1->ignore_best = dom;
1201  }
1202  p1 = p1->next[ord];
1203  }
1204 
1205  /* If there is any point dominated in every dimension, it is removed and
1206  the number of total points is updated. */
1207  if (count[dim] > 0) {
1208  p1 = list->prev[0];
1209  while (p1->x) {
1210  if (p1->ignore == dim) {
1211  delete_dominated(p1, dim);
1212  c--;
1213  }
1214  p1 = p1->prev[0];
1215  }
1216  }
1217  return c;
1218 }
1219 
1220 
1221 /*
1222  Tries to find a good order to process the objectives.
1223 
1224  This algorithm tries to maximize the number of dominated points
1225  dominated in more dimensions. For example, for a problem with d
1226  dimensions, an order with 20 points dominated from dimension 0 to
1227  dimension d-1 is prefered to an order of objectives in which the
1228  number of points dominated from dimension 0 to d-1 is 10. An order
1229  with the same number of points dominated up to dimension d-1 as a
1230  second order is prefered if it has more points dominated up to
1231  dimension d-2 than the second order. */
1232 static int define_order(dlnode_t *list, int dim, int c, int *order)
1233 {
1234  dlnode_t *p;
1235 
1236  // order - keeps the current order of objectives
1237 
1238  /* best_order - keeps the current best order for the
1239  objectives. At the end, this array (and the array order) will
1240  have the best order found, to process the objectives.
1241 
1242  This array keeps the indexes of the objectives, where
1243  best_order[0] keeps the index of the first objective,
1244  best_order[1] keeps the index of the second objective and so on. */
1245  int *best_order = malloc(dim * sizeof(int));
1246 
1247  /* count - keeps the counting of the dominated points
1248  corresponding to the order of objectives in 'order'.
1249 
1250  When it's found that a point is dominated at most, for the
1251  first four dimensions, then count[3] is incremented. So,
1252  count[i] is incremented every time it's found a point that is
1253  dominated from dimension 0 to i, but not in dimension i+1. */
1254  int *count = malloc(dim * sizeof(int));
1255 
1256  /* keeps the best counting of the dominated points (that is
1257  obtained using the order in best_order). */
1258  int *best_count = malloc(dim * sizeof(int));
1259 
1260  int i, j, k;
1261 
1262  for (i = 0; i < dim; i++) {
1263  best_order[i] = order[i] = i;
1264  best_count[i] = count[i] = 0;
1265  }
1266 
1267  // determines the number of dominated points in the original order.
1268  // c - total number of points excluding points totally dominated
1269  c = determine_ndomf(list, dim-1, c, order, count);
1270 
1271  /* the best order so far is the original order, so it's necessary
1272  to register the number of points dominated in the best
1273  order. */
1274  for (i = 0; i < dim; i++) {
1275  best_count[i] = count[i];
1276  }
1277 
1278  /* Objectives are chosen from highest to lowest. So we start
1279  defining which is the objective in position dim-1 and then
1280  which is the objective in position dim, and so on. The
1281  objective chosen to be in position i is chosen in a way to
1282  maximize the number of dominated points from dimension 0 to
1283  i-1. So, this cycle, selects a position i, and then we find
1284  the objective (from the remaining objectives that haven't a
1285  position yet, the objectives that are in positions lower or
1286  equal to i) that by being in position i maximizes the number of
1287  points dominated from dimension 0 to i-1. */
1288  for (i = dim - 1; i > 2; i--) {
1289  /* This cycle, in every iteration, assigns a different
1290  objective to position i. It's important to notice
1291  that if we want to maximize the number of dominated
1292  points from dimension 0 to i-1, when we want to now if
1293  an objective k in position i is the one that maximizes
1294  it, it doesn't matter the order of the objectives in
1295  positions lower than i, the number of dominated points
1296  from dimension 0 to i-1 will always be the same, so
1297  it's not necessary to worry about the order of those
1298  objectives.
1299 
1300  When this cycle starts, 'order' has the original order and
1301  so 'count' has the number of points dominated from 0
1302  to every k, where k < dim or 'order' has the last
1303  order of objectives used to calculate the best
1304  objective to put in position i+1 that maximizes the
1305  number of dominated points from dimension 0 to i and
1306  so 'count' has the number of points dominated from
1307  dimension 0 to every k, where k < dim, that was
1308  calculated previously.
1309 
1310  There on, it is not necessary to calculate the number of
1311  dominated points from dimension 0 to i-1 with the actual
1312  objective in position i (order[i]), because this value was
1313  previously calculated and so it is only necessary to
1314  calculate the number of dominated points when the current
1315  objectives in order[k], where k < i, are in position i. */
1316  for (j = 0; j < i; j++) {
1317  int aux = order[i];
1318  order[i] = order[j];
1319  order[j] = aux;
1320 
1321  /* Determine the number of dominated points from dimension
1322  0 to k, where k < i (the number of points dominated
1323  from dimension 0 to t, where t >= i, is already known
1324  from previous calculations) with a different objective
1325  in position i. */
1326  determine_ndom(list, i-1, order, count);
1327 
1328  /* If the order in 'order' is better than the previously
1329  best order, than the actual order is now the best. An
1330  order is better than another if the number of dominated
1331  points from dimension 0 to i-1 is higher. If this
1332  number is equal, then the best is the one that has the
1333  most dominated points from dimension 0 to i-2. If this
1334  number is equal, than the last order considered the
1335  best, still remains the best order so far. */
1336  if (best_count[i-1] < count[i-1]
1337  || (best_count[i-1] == count[i-1]
1338  && best_count[i-2] < count[i-2])) {
1339  for (k = 0; k <= i; k++) {
1340  best_count[k] = count[k];
1341  best_order[k] = order[k];
1342  }
1343  p = list->prev[0];
1344  while (p != list) {
1345  p->ignore_best = p->ignore;
1346  p = p->prev[0];
1347  }
1348  }
1349  }
1350 
1351  /*
1352  If necessary, update 'order' with the best order so far and
1353  the corresponding number of dominated points. In this way,
1354  in the next iteration it is not necessary to recalculate the
1355  number of dominated points from dimension 0 to i-2, when in
1356  position i-1 is the objective that is currently in position
1357  i-1, in the best order so far (best_order[i-1]).
1358  */
1359  if (order[i] != best_order[i]) {
1360  for (j = 0; j <= i; j++) {
1361  count[j] = best_count[j];
1362  order[j] = best_order[j];
1363  }
1364  p = list->prev[0];
1365  /*
1366  The information about a point being dominated is updated
1367  because, this way, in some cases it is not necessary to
1368  find out (again) if a point is dominated.
1369  */
1370  while (p != list) {
1371  p->ignore = p->ignore_best;
1372  p = p->prev[0];
1373  }
1374  }
1375 
1376  }
1377 
1378  free(count);
1379  free(best_count);
1380  free(best_order);
1381  return c;
1382 }
1383 
1384 /*
1385  Reorders the reference point's objectives according to an order 'order'.
1386 */
1387 static void reorder_reference(double *reference, int d, int *order)
1388 {
1389  int j;
1390  double *tmp = (double *) malloc(d * sizeof(double));
1391  for (j = 0; j < d; j++) {
1392  tmp[j] = reference[j];
1393  }
1394  for (j = 0; j < d; j++) {
1395  reference[j] = tmp[order[j]];
1396  }
1397  free(tmp);
1398 }
1399 
1400 /*
1401  Reorders the dimensions for every point according to an order.
1402 */
1403 void reorder_list(dlnode_t *list, int d, int *order)
1404 {
1405  int j;
1406  double *x;
1407  double *tmp = (double *) malloc(d * sizeof(double));
1408  dlnode_t **prev = (dlnode_t **) malloc(d * sizeof(dlnode_t *));
1409  dlnode_t **next = (dlnode_t **) malloc(d * sizeof(dlnode_t *));
1410  dlnode_t *p;
1411 
1412  for(j = 0; j < d; j++) {
1413  prev[j] = list->prev[j];
1414  next[j] = list->next[j];
1415  }
1416 
1417  for(j = 0; j < d; j++) {
1418  list->prev[j] = prev[order[j]];
1419  list->next[j] = next[order[j]];
1420  }
1421  p = list->next[0];
1422 
1423  while (p != list) {
1424  p->ignore = 0;
1425  x = p->x;
1426  for(j = 0; j < d; j++) {
1427  tmp[j] = x[j];
1428  prev[j] = p->prev[j];
1429  next[j] = p->next[j];
1430  }
1431 
1432  for(j = 0; j < d; j++) {
1433  x[j] = tmp[order[j]];
1434  p->prev[j] = prev[order[j]];
1435  p->next[j] = next[order[j]];
1436  }
1437  p = p->next[0];
1438  }
1439  free(tmp);
1440  free(prev);
1441  free(next);
1442 }
1443 #endif
1444 
1445 double fpli_hv(double *data, int d, int n, const double *ref)
1446 {
1447  dlnode_t *list;
1448  double hyperv;
1449  double * bound = NULL;
1450 
1451 #if VARIANT >= 3
1452  int i;
1453 
1454  bound = malloc (d * sizeof(double));
1455  for (i = 0; i < d; i++) bound[i] = -DBL_MAX;
1456 #endif
1457 
1458  tree = avl_alloc_tree ((avl_compare_t) compare_tree_asc,
1459  (avl_freeitem_t) NULL);
1460 
1461  list = setup_cdllist(data, d, n);
1462 
1463  n = filter(list, d, n, ref);
1464  if (n == 0) {
1465  hyperv = 0.0;
1466  } else if (n == 1) {
1467  dlnode_t * p = list->next[0];
1468  hyperv = 1;
1469  for (i = 0; i < d; i++)
1470  hyperv *= ref[i] - p->x[i];
1471  } else {
1472  hyperv = hv_recursive(list, d-1, n, ref, bound);
1473  }
1474  /* Clean up. */
1475  free_cdllist (list);
1476  free (tree); /* The nodes are freed by free_cdllist (). */
1477  free (bound);
1478 
1479  return hyperv;
1480 }
1481 
1482 #ifdef EXPERIMENTAL
1483 
1484 #include "timer.h" /* FIXME: Avoid calling Timer functions here. */
1485 double fpli_hv_order(double *data, int d, int n, const double *ref, int *order,
1486  double *order_time, double *hv_time)
1487 {
1488  dlnode_t *list;
1489  double hyperv;
1490  double * bound = NULL;
1491  double * ref_ord = (double *) malloc(d * sizeof(double));
1492 
1493 #if VARIANT >= 3
1494  int i;
1495 
1496  bound = malloc (d * sizeof(double));
1497  for (i = 0; i < d; i++) bound[i] = -DBL_MAX;
1498 #endif
1499 
1500  tree = avl_alloc_tree ((avl_compare_t) compare_tree_asc,
1501  (avl_freeitem_t) NULL);
1502 
1503  list = setup_cdllist(data, d, n);
1504 
1505  if (d > 3) {
1506  n = define_order(list, d, n, order);
1507  reorder_list(list, d, order);
1508 
1509  // copy reference so it will be unchanged for the next data sets.
1510  for (i = 0; i < d; i++)
1511  ref_ord[i] = ref[i];
1512 
1513  reorder_reference(ref_ord, d, order);
1514 
1515  } else {
1516  for(i = 0; i < d; i++)
1517  ref_ord[i] = ref[i];
1518  }
1519 
1520  *order_time = Timer_elapsed_virtual ();
1521  Timer_start();
1522 
1523  n = filter(list, d, n, ref_ord);
1524  if (n == 0) {
1525  hyperv = 0.0;
1526  } else if (n == 1) {
1527  hyperv = 1;
1528  dlnode_t * p = list->next[0];
1529  for (i = 0; i < d; i++)
1530  hyperv *= ref[i] - p->x[i];
1531  } else {
1532  hyperv = hv_recursive(list, d-1, n, ref, bound);
1533  }
1534  /* Clean up. */
1535  free_cdllist (list);
1536  free (tree); /* The nodes are freed by free_cdllist (). */
1537  free (bound);
1538  free (ref_ord);
1539 
1540  *hv_time = Timer_elapsed_virtual ();
1541 
1542  return hyperv;
1543 }
1544 #endif
Generated by   doxygen 1.8.9.1