diff --git a/c/Makefile b/c/Makefile index 13dfd403f..67733f9ed 100644 --- a/c/Makefile +++ b/c/Makefile @@ -97,6 +97,7 @@ HDRS = avl.h \ hv_priv.h \ hv3d_priv.h \ hv4d_priv.h \ + hvc4d_priv.h \ igd.h \ io.h \ io_priv.h \ @@ -204,7 +205,7 @@ eaf.o: eaf.h io.h bit_array.h cvector.h eaf3d.o: eaf.h io.h bit_array.h cvector.h avl.h eaf_main.o: cmdline.h io.h eaf.h bit_array.h cvector.h epsilon.o: cmdline.h io.h epsilon.h nondominated.h sort.h avl_tiny.h -hv.o: hv.h hv_priv.h hv4d_priv.h sort.h libmoocore-config.h +hv.o: hv.h hv_priv.h hvc4d_priv.h sort.h libmoocore-config.h hv3dplus.o: hv_priv.h sort.h hv3d_priv.h avl_tiny.h hv4d.o: hv4d_priv.h hv_priv.h sort.h hv_contrib.o: hv.h libmoocore-config.h nondominated.h sort.h avl_tiny.h diff --git a/c/NEWS.md b/c/NEWS.md index 03d7f251d..36ee067d7 100644 --- a/c/NEWS.md +++ b/c/NEWS.md @@ -12,6 +12,8 @@ * Reorganize igd.h so that helper functions are inlined and more loops are vectorized. * Change the type of `minmax` from `signed char *` to `int *` to help vectorization. * Fix bug in `epsilon_mult()` with mixed min-max objectives. + * `fpli_hv()` calculates 4D contributions as the base case, which is + significantly faster for more than four dimensions. (Andreia P. Guerreiro) ## 0.17.0 diff --git a/c/hv.c b/c/hv.c index ac043e285..b02336523 100644 --- a/c/hv.c +++ b/c/hv.c @@ -35,7 +35,7 @@ #include "common.h" #include "hv.h" #define HV_RECURSIVE -#include "hv4d_priv.h" +#include "hvc4d_priv.h" #define STOP_DIMENSION 3 // stop on dimension 4. #define MAX_ROWS_HV_INEX 15 @@ -64,7 +64,11 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, dimension_t d_stop = d - STOP_DIMENSION; size_t n = *size; - dlnode_t * head = malloc((n+1) * sizeof(*head)); + // Reserve space for main CDLL used by hv_recursive() and for the auxiliary + // list used by onec4dplusU(). The main CDLL will store all points + 1 + // sentinel. The auxiliary list will store at most n - 1 points + 1 + // sentinel. + dlnode_t * head = malloc((n + 1 + n) * sizeof(*head)); size_t i = 1; for (size_t j = 0; j < n; j++) { /* Filters those points that do not strictly dominate the reference @@ -85,9 +89,10 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, // We need space in r_next and r_prev for dimension 5 and above (d_stop - 1). head->r_next = malloc(2 * (d_stop - 1) * (n+1) * sizeof(head)); head->r_prev = head->r_next + (d_stop - 1) * (n+1); - // We only need space in area and vol for dimension 4 and above. - head->area = malloc(2 * d_stop * (n+1) * sizeof(*data)); - head->vol = head->area + d_stop * (n+1); + // We only need space in area and vol for dimension 4 and above (d-stop). + // Also reserve space for n-1 auxiliary 3D points used by onec4dplusU(). + head->area = malloc((2 * d_stop * n + 3 * (n-1)) * sizeof(*data)); + head->vol = head->area + d_stop * n; head->x = NULL; // head contains no data head->ignore = 0; // should never get used @@ -96,7 +101,13 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, // Link head and list4d; head is not used by HV4D, so next[0] and prev[0] // should remain untouched. head->next[0] = list4d; - head->prev[0] = list4d; // Save it twice so we can use assert() later. + + // Setup the auxiliary list used in onec4dplusU(). + dlnode_t * list_aux = head + n + 1; + // Setup the auxiliary 3D points used in onec4dplusU(). + list_aux->vol = head->vol + d_stop * n; + head->prev[0] = list_aux; + list_aux->next[0] = list4d; for (i = 1; i <= n; i++) { // Shift x because qsort() cannot take the dimension to sort as an argument. @@ -104,33 +115,24 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, head[i].ignore = 0; head[i].r_next = head->r_next + i * (d_stop - 1); head[i].r_prev = head->r_prev + i * (d_stop - 1); - head[i].area = head->area + i * d_stop; - head[i].vol = head->vol + i * d_stop; + head[i].area = head->area + (i - 1) * d_stop; + head[i].vol = head->vol + (i - 1) * d_stop; } + // Make sure they are not used. + head->area = NULL; + head->vol = NULL; dlnode_t ** scratch = malloc(n * sizeof(*scratch)); for (i = 0; i < n; i++) scratch[i] = head + 1 + i; - int j = d_stop - 2; - while (true) { + for (int j = d_stop - 2; j >= 0; j--) { /* FIXME: replace qsort() by something better: https://github.com/numpy/x86-simd-sort https://github.com/google/highway/tree/52a2d98d07852c5d69284e175666e5f8cc7d8285/hwy/contrib/sort */ // Sort each dimension independently. qsort(scratch, n, sizeof(*scratch), compare_node); - if (j == -1) { - (list4d+1)->next[1] = scratch[0]; - scratch[0]->prev[1] = list4d+1; - for (i = 1; i < n; i++) { - scratch[i-1]->next[1] = scratch[i]; - scratch[i]->prev[1] = scratch[i-1]; - } - scratch[n-1]->next[1] = list4d+2; - (list4d+2)->prev[1] = scratch[n-1]; - break; - } head->r_next[j] = scratch[0]; scratch[0]->r_prev[j] = head; for (i = 1; i < n; i++) { @@ -139,21 +141,34 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } scratch[n-1]->r_next[j] = head; head->r_prev[j] = scratch[n-1]; - j--; // Consider next objective (in reverse order). for (i = 1; i <= n; i++) head[i].x--; } - // Reset x to point to the first objective. - for (i = 1; i <= n; i++) - head[i].x -= STOP_DIMENSION; + for (int j = 1; j >= 0; j--) { + // Sort each dimension independently. + qsort(scratch, n, sizeof(*scratch), compare_node); + (list4d+1)->next[j] = scratch[0]; + scratch[0]->prev[j] = list4d+1; + 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] = list4d+2; + (list4d+2)->prev[j] = scratch[n-1]; + if (j > 0) { + // Consider next objective (in reverse order). + for (i = 1; i <= n; i++) + head[i].x--; + } else { + // Reset x to point to the first objective. + for (i = 1; i <= n; i++) + head[i].x -= STOP_DIMENSION - 1; + } + } free(scratch); - // Make sure it is not used. - ASAN_POISON_MEMORY_REGION(head->area, sizeof(*data) * d_stop); - ASAN_POISON_MEMORY_REGION(head->vol, sizeof(*data) * d_stop); - finish: *size = n; return head; @@ -161,11 +176,11 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, static void fpli_free_cdllist(dlnode_t * head) { - assert(head->next[0] == head->prev[0]); + assert(head->next[0] == head->prev[0]->next[0]); free_cdllist(head->next[0]); // free 4D sentinels free(head->r_next); - free(head->area); - free(head); + free(head[1].area); // Free ->area, ->vol and list_aux->vol (x_aux used by 4D basecase). + free(head); // Free main CDLL and list_aux (4D basecase). } static inline void @@ -180,6 +195,34 @@ update_bound(double * restrict bound, const double * restrict x, dimension_t dim } } +static void +delete_4d(dlnode_t * restrict nodep) +{ + nodep->prev[1]->next[1] = nodep->next[1]; + nodep->next[1]->prev[1] = nodep->prev[1]; +} + +static void +delete_3d(dlnode_t * restrict nodep) +{ + nodep->prev[0]->next[0] = nodep->next[0]; + nodep->next[0]->prev[0] = nodep->prev[0]; +} + +static void +reinsert_4d(dlnode_t * restrict nodep) +{ + nodep->prev[1]->next[1] = nodep; + nodep->next[1]->prev[1] = nodep; +} + +static void +reinsert_3d(dlnode_t * restrict nodep) +{ + nodep->prev[0]->next[0] = nodep; + nodep->next[0]->prev[0] = nodep; +} + static void delete_dom(dlnode_t * restrict nodep, dimension_t dim) { @@ -190,9 +233,8 @@ delete_dom(dlnode_t * restrict nodep, dimension_t dim) nodep->r_prev[d]->r_next[d] = nodep->r_next[d]; nodep->r_next[d]->r_prev[d] = nodep->r_prev[d]; } - // Dimension 4. - nodep->prev[1]->next[1] = nodep->next[1]; - nodep->next[1]->prev[1] = nodep->prev[1]; + delete_4d(nodep); + delete_3d(nodep); } static void @@ -213,9 +255,8 @@ reinsert_nobound(dlnode_t * restrict nodep, dimension_t dim) nodep->r_prev[d]->r_next[d] = nodep; nodep->r_next[d]->r_prev[d] = nodep; } - // Dimension 4. - nodep->prev[1]->next[1] = nodep; - nodep->next[1]->prev[1] = nodep; + reinsert_4d(nodep); + reinsert_3d(nodep); } static void @@ -226,15 +267,14 @@ reinsert(dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) } static double -fpli_hv4d(dlnode_t * restrict list, size_t c _attr_maybe_unused) +fpli_onec4d(dlnode_t * restrict list, size_t c _attr_maybe_unused, dlnode_t * restrict the_point) { ASSUME(c > 1); - assert(list->next[0] == list->prev[0]); + assert(list->next[0] == list->prev[0]->next[0]); dlnode_t * restrict list4d = list->next[0]; - // hv4dplusU() will change the sentinels for 3D, so we need to reset them. - reset_sentinels_3d(list4d); - double hv = hv4dplusU(list4d); - return hv; + dlnode_t * restrict list_aux = list->prev[0]; + double contrib = onec4dplusU(list4d, list_aux, the_point); + return contrib; } _attr_optimize_finite_and_associative_math // Required for auto-vectorization: https://gcc.gnu.org/PR122687 @@ -248,14 +288,6 @@ one_point_hv(const double * restrict x, const double * restrict ref, dimension_t return hv; } -_attr_optimize_finite_and_associative_math -static inline void -upper_bound(double * restrict dest, const double * restrict a, const double * restrict b, dimension_t dim) -{ - for (dimension_t i = 0; i < dim; i++) - dest[i] = MAX(a[i], b[i]); -} - _attr_optimize_finite_and_associative_math static double hv_two_points(const double * restrict x1, const double * restrict x2, @@ -360,14 +392,6 @@ hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, const double * restrict ref, double * restrict bound) { ASSUME(c > 1); - ASSUME(dim >= STOP_DIMENSION); - if (dim == STOP_DIMENSION) { - /*--------------------------------------- - base case of dimension 4 - --------------------------------------*/ - return fpli_hv4d(list, c); - } - ASSUME(dim > STOP_DIMENSION); /* ------------------------------------------------------ General case for dimensions higher than 4D ------------------------------------------------------ */ @@ -434,7 +458,15 @@ hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, DEBUG1(debug_counter[1]++); hypera = p1_prev->area[d_stop]; } else { - hypera = hv_recursive(list, dim - 1, c, ref, bound); + ASSUME(dim - 1 >= STOP_DIMENSION); + if (dim - 1 == STOP_DIMENSION) { + // base case of dimension 4. + hypera = fpli_onec4d(list, c, p1); + // hypera only has the contribution of p1. + hypera += p1_prev->area[d_stop]; + } else { + hypera = hv_recursive(list, dim - 1, c, ref, bound); + } /* At this point, p1 is the point with the highest value in dimension dim in the list: If it is dominated in dimension dim-1, then it is also dominated in dimension dim. */ @@ -503,7 +535,16 @@ fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, p1->vol[d_stop] = hyperv; assert(p1->ignore == 0); c++; - double hypera = hv_recursive(list, dim - 1, c, ref, bound); + double hypera; + ASSUME(dim - 1 >= STOP_DIMENSION); + if (dim - 1 == STOP_DIMENSION) { + // base case of dimension 4. + hypera = fpli_onec4d(list, c, p1); + // hypera only has the contribution of p1. + hypera += p1_prev->area[d_stop]; + } else { + hypera = hv_recursive(list, dim - 1, c, ref, bound); + } /* At this point, p1 is the point with the highest value in dimension dim in the list: If it is dominated in dimension dim-1, then it is also dominated in dimension dim. */ diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index f92b8a340..7c51c6c0b 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -101,7 +101,11 @@ static bool restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict newp) { // FIXME: This is the most expensive function in the HV4D+ algorithm. +#ifdef HV_RECURSIVE + const double newx[] = { newp->x[0], newp->x[1], newp->x[2] }; +#else const double newx[] = { newp->x[0], newp->x[1], newp->x[2], newp->x[3] }; +#endif assert(list+1 == list->next[0]); dlnode_t * closest0 = list+1; dlnode_t * closest1 = list; @@ -126,7 +130,12 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n // if (weakly_dominates(px, newx, 3)) { // Slower if (p_leq_new_0 & p_leq_new_1 & p_leq_new_2) { //new->ndomr++; +#ifdef HV_RECURSIVE + // When called by onec4dplusU, px may be just a 3-dimensional vector. + assert(weakly_dominates(px, newx, 3)); +#else assert(weakly_dominates(px, newx, 4)); +#endif return false; } @@ -173,7 +182,11 @@ one_contribution_3d(dlnode_t * restrict newp) double volume = 0; double lastz = newx[2]; dlnode_t * p = newp->next[0]; +#if !defined(HV_RECURSIVE) && HV_DIMENSION == 4 assert(!weakly_dominates(p->x, newx, 4)); +#else + assert(!weakly_dominates(p->x, newx, 3)); +#endif while (true) { const double * px = p->x; volume += area * (px[2] - lastz); @@ -219,7 +232,7 @@ one_contribution_3d(dlnode_t * restrict newp) /* Compute the hypervolume indicator in d=4 by iteratively computing the one contribution problem in d=3. */ -static double +static inline double hv4dplusU(dlnode_t * list) { assert(list+2 == list->prev[0]); diff --git a/c/hv_priv.h b/c/hv_priv.h index 352c15c58..7ee2d1ee5 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -325,4 +325,14 @@ compute_area_no_inners(const double * px, const dlnode_t * q, uint_fast8_t i) ASSUME(i == 0 || i == 1); return compute_area_simple(px, q, q->cnext[i], i); } + +_attr_optimize_finite_and_associative_math +static inline void +upper_bound(double * restrict dest, const double * restrict a, + const double * restrict b, dimension_t dim) +{ + for (dimension_t i = 0; i < dim; i++) + dest[i] = MAX(a[i], b[i]); +} + #endif // _HV_PRIV_H diff --git a/c/hvc4d_priv.h b/c/hvc4d_priv.h new file mode 100644 index 000000000..08f5b793e --- /dev/null +++ b/c/hvc4d_priv.h @@ -0,0 +1,317 @@ +#ifndef _HVC4D_PRIV_H +#define _HVC4D_PRIV_H +/****************************************************************************** + Compute hypervolume contributions in 4D. Currently only used by hv.c + ------------------------------------------------------------------------------ + + Copyright (C) 2025--2026 + Andreia P. Guerreiro + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + ------------------------------------------------------------------------------ + + Reference: + + [1] Andreia P. Guerreiro and Carlos M. Fonseca. Computing and Updating + Hypervolume Contributions in Up to Four Dimensions. IEEE Transactions on + Evolutionary Computation, 22(3):449–463, June 2018. + + ------------------------------------------------------------------------------ + +******************************************************************************/ + +#include "hv4d_priv.h" + +typedef struct lex_sort_buffer { + const double ** scratch; + size_t capacity; +} lex_sort_buffer_t; + +static inline int +lex_sort_equal_z_and_setup_nodes(lex_sort_buffer_t * restrict buff, + dlnode_t * restrict newp_aux, + const double * restrict x_aux, size_t n) +{ + const double ** scratch = buff->scratch; + if (buff->capacity < n) { + scratch = realloc(scratch, n * sizeof(*scratch)); + if (!scratch) { + buff->scratch = NULL; + return -1; + } + buff->scratch = scratch; + buff->capacity = n; + } + + const double * x = x_aux; + for (size_t i = 0; i < n; i++){ + scratch[i] = x; + x += 3; + } + + qsort(scratch, n, sizeof(*scratch), cmp_ppdouble_asc_rev_2d); + + for (size_t i = 0; i < n; i++) { + newp_aux->x = scratch[i]; + assert(i == 0 || lexicographic_less_2d(scratch[i-1], scratch[i])); + newp_aux++; + } + return 0; +} + + +/** + Assumes that the HV3D+ data structure is reconstructed up to z <= newp->x[2]. + Sweeps the points in the data structure in ascending order of y-coordinate, includes newp + in the z list, sets up the "closest" data of newp and, if equalz=true (i.e., all points in + the data structure have x[2] == newp->x[2]), it also updates the "closest" data of the point + that follows newp according to the lexicographic order (z,y,x). +*/ +static inline bool +continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict newp, bool equalz) +{ + const double newx[] = { newp->x[0], newp->x[1], newp->x[2] }; + assert(list+1 == list->next[0]); + dlnode_t * p = (list+1)->cnext[1]; + + // find newp->closest[0] + while (p->x[1] < newx[1]){ + assert(p->x[2] <= newx[2]); + p = p->cnext[1]; + } + + if (p->x[1] != newx[1] || p->x[0] > newx[0]) + p = p->cnext[0]; + + assert(lexicographic_less_2d(p->x, newx)); + assert(lexicographic_less_2d(newx, p->cnext[1]->x)); + + // Check if newp is dominated. + if (weakly_dominates(p->x, newx, 3)) { + return false; + } + dlnode_t * lex_prev = p; // newp->closest[0] = lex_prev + p = lex_prev->cnext[1]; + + // if all points in the list have z-coordinate equal to px[2] + if (equalz) { + assert(p->cnext[0]->x[1] < newx[1]); + assert(p->x[1] >= newx[1]); + // Check if newp dominates points in the list. + while (p->x[0] >= newx[0]) { + remove_from_z(p); + p = p->cnext[1]; + } + + assert(p->x[0] < newx[0]); + // update the closest of points in the list (max. 1), if needed + if (p != list) + p->closest[0] = newp; + + //setup z list + assert(p == list || lex_prev->next[0] == p); + newp->prev[0] = lex_prev; + newp->next[0] = lex_prev->next[0]; + + newp->closest[1] = list; + + } else { + //setup z list + newp->prev[0] = list->prev[0]->prev[0]; + newp->next[0] = list->prev[0]; + + // find newp->closest[1] + while (p->x[0] >= newx[0]){ + assert(p->x[2] <= newx[2]); + assert(!weakly_dominates(newx, p->x, 3)); + p = p->cnext[1]; + } + newp->closest[1] = p; + } + + newp->closest[0] = lex_prev; + // update cnext + lex_prev->cnext[1] = newp; + newp->cnext[0] = lex_prev; + newp->cnext[1] = p; + p->cnext[0] = newp; + + //update z list + newp->next[0]->prev[0] = newp; + newp->prev[0]->next[0] = newp; + + return true; +} + +/** + Compute the hv contribution of "the_point" in d=4 by iteratively computing + the one contribution problem in d=3. */ +static inline double +onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, + dlnode_t * restrict the_point) +{ + // FIXME: This is never triggered. + assert(the_point->ignore < 3); + if (the_point->ignore >= 3) { + return 0; + } + + assert(list+2 == list->prev[0]); + assert(list+2 == list->prev[1]); + assert(list+1 == list->next[1]); + + const dlnode_t * const last = list+2; + dlnode_t * const z_first = (list+1)->next[0]; + dlnode_t * const z_last = last->prev[0]; + + double * x_aux = list_aux->vol; + dlnode_t * newp_aux = list_aux+1; // list_aux is a sentinel + + reset_sentinels_3d(list); + restart_list_y(list); + + const double * the_point_x = the_point->x; + // Setup the 3D base only if there are any points leq than the_point_x[3]. + if ((list+1)->next[1] != the_point || the_point_x[3] == the_point->next[1]->x[3]) { + bool done_once = false; + // Set the_point->ignore=3 so the loop will skip it, but restore its + // value after the loop. + dimension_t the_point_ignore = the_point->ignore; + the_point->ignore = 3; + dlnode_t * newp = z_first; + assert(newp != last); + + // PART 1: Setup 2D base of the 3D base + while (newp->x[2] <= the_point_x[2]) { + const double * newpx = newp->x; + if (newpx[3] <= the_point_x[3] && newp->ignore < 3) { + if (weakly_dominates(newpx, the_point_x, 3)) { + assert(the_point->ignore == 3); + // Restore modified links (z list). + (list+1)->next[0] = z_first; + (list+2)->prev[0] = z_last; + return 0; + } + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + upper_bound(x_aux, newpx, the_point_x, 3); + newp_aux->x = x_aux; + x_aux += 3; + + if (continue_base_update_z_closest(list, newp_aux, done_once)) { + newp_aux++; + done_once = true; + } + } + newp = newp->next[0]; + } + the_point->ignore = the_point_ignore; + + // PART 2: Setup the remainder of the 3D base + int c = 0; + // FIXME: Could we keep this buffer in list_aux? Or pass it down as an argument? + lex_sort_buffer_t sort_buffer = { .scratch = NULL, .capacity = 0 }; + while (newp != last) { + const double * newpx = newp->x; + if (newpx[3] <= the_point_x[3] && newp->ignore < 3) { + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + upper_bound(x_aux, newpx, the_point_x, 3); + x_aux += 3; + c++; + } + + if (c > 0 && newp->next[0]->x[2] > newpx[2]) { + if (c == 1) { + newp_aux->x = x_aux - 3; + continue_base_update_z_closest(list, newp_aux, false); + newp_aux++; + } else { + // all points with equal z-coordinate will be added to the data structure in lex order + lex_sort_equal_z_and_setup_nodes(&sort_buffer, newp_aux, x_aux - 3*c, c); + continue_base_update_z_closest(list, newp_aux, false); + const double * prevp_aux_x = newp_aux->x; + newp_aux++; + + for (int i = 1; i < c; i++) { + assert(newpx[2] == newp_aux->x[2]); // all c points have equal z + assert(prevp_aux_x[1] <= newp_aux->x[1]); // due to lexsort + if (newp_aux->x[0] < prevp_aux_x[0]){ + // if newp_aux is not dominated by prevp + continue_base_update_z_closest(list, newp_aux, false); + } + prevp_aux_x = newp_aux->x; + newp_aux++; + } + } + c = 0; + } + newp = newp->next[0]; + } + free(sort_buffer.scratch); + } + + dlnode_t * newp = the_point->next[1]; + while (newp->x[3] <= the_point_x[3]) + newp = newp->next[1]; + + dlnode_t * tp_prev_z = the_point->prev[0]; + dlnode_t * tp_next_z = the_point->next[0]; + // This call should always return true. + _attr_maybe_unused bool restart_base_setup_z_and_closest_result = + restart_base_setup_z_and_closest(list, the_point); + assert(restart_base_setup_z_and_closest_result); + double volume = one_contribution_3d(the_point); + the_point->prev[0] = tp_prev_z; + the_point->next[0] = tp_next_z; + + assert(volume > 0); + double height = newp->x[3] - the_point_x[3]; + // It cannot be zero because we exited the loop above. + assert(height > 0); + double hv = volume * height; + + // PART 3: Update the 3D contribution. + while (newp != last) { + const double * newpx = newp->x; + if (weakly_dominates(newpx, the_point_x, 3)) + break; + if (newp->ignore < 3) { + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + upper_bound(x_aux, newpx, the_point_x, 3); + newp_aux->x = x_aux; + x_aux += 3; + if (restart_base_setup_z_and_closest(list, newp_aux)) { + // newp was not dominated by something else. + double newp_v = one_contribution_3d(newp_aux); + assert(newp_v > 0); + volume -= newp_v; + + add_to_z(newp_aux); + update_links(list, newp_aux); + newp_aux++; + } + + if (weakly_dominates(the_point_x, newpx, 3)){ + newp->ignore = 3; + } + } + + // FIXME: If newp was dominated, can we accumulate the height and update + // hv later? AG: Yes + height = newp->next[1]->x[3] - newpx[3]; + assert(height >= 0); + hv += volume * height; + + newp = newp->next[1]; + } + + // Restore z list + (list+1)->next[0] = z_first; + (list+2)->prev[0] = z_last; + + return hv; +} + +#endif // _HVC4D_PRIV_H diff --git a/c/libhv.mk b/c/libhv.mk index faa368ab9..53f34492b 100644 --- a/c/libhv.mk +++ b/c/libhv.mk @@ -1,6 +1,6 @@ # -*- Makefile-gmake -*- LIBHV_SRCS = hv.c hv3dplus.c hv4d.c hvc3d.c hv_contrib.c -LIBHV_HDRS = hv.h hv_priv.h hv3d_priv.h hv4d_priv.h libmoocore-config.h +LIBHV_HDRS = hv.h hv_priv.h hv3d_priv.h hv4d_priv.h hvc4d_priv.h libmoocore-config.h LIBHV_OBJS = $(LIBHV_SRCS:.c=.o) HV_LIB = fpli_hv.a diff --git a/c/sort.h b/c/sort.h index 3f78fd5a7..9a59b2843 100644 --- a/c/sort.h +++ b/c/sort.h @@ -54,6 +54,12 @@ lexicographic_less_3d(const double * restrict a, const double * restrict b) return a[2] < b[2] || (a[2] == b[2] && (a[1] < b[1] || (a[1] == b[1] && a[0] <= b[0]))); } +static inline bool +lexicographic_less_2d(const double * restrict a, const double * restrict b) +{ + return a[1] < b[1] || (a[1] == b[1] && a[0] <= b[0]); +} + static inline bool all_equal_double(const double * restrict a, const double * restrict b, dimension_t dim) { @@ -110,9 +116,8 @@ cmp_pdouble_asc_rev_2d(const void * restrict pa, const void * restrict pb) { const double * restrict a = (const double *) pa; const double * restrict b = (const double *) pb; - const double ax = a[0], ay = a[1], bx = b[0], by = b[1]; - int cmpx = cmp_double_asc(ax, bx); - int cmpy = cmp_double_asc(ay, by); + int cmpx = cmp_double_asc(a[0], b[0]); + int cmpy = cmp_double_asc(a[1], b[1]); return cmpy ? cmpy : cmpx; } static inline int diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index 1c672b03c..5bd2b1346 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -465,14 +465,15 @@ def hypervolume( :math:`O(m 2^{n})` time and :math:`O(n\cdot m)` space complexity, where :math:`m` is the dimension of the points, but it is very fast for such small sets. For larger number of points, it uses a recursive algorithm - :footcite:p:`FonPaqLop06:hypervolume` with HV4D\ :sup:`+` as the base case, - resulting in a :math:`O(n^{m-2})` time complexity and :math:`O(n)` space - complexity in the worst-case. Experimental results show that the pruning - techniques used may reduce the time complexity even further. The original - proposal :footcite:p:`FonPaqLop06:hypervolume` had the HV3D algorithm as - the base case, giving a time complexity of :math:`O(n^{m-2} \log n)`. - Andreia P. Guerreiro enhanced the numerical stability of the recursive - algorithm by avoiding floating-point comparisons of partial hypervolumes. + :footcite:p:`FonPaqLop06:hypervolume` that computes 4D contributions + :footcite:p:`GueFon2017hv4d` as the base case, resulting in a + :math:`O(n^{m-2})` time complexity and :math:`O(n)` space complexity in the + worst-case. Experimental results show that the pruning techniques used may + reduce the time complexity even further. The original proposal + :footcite:p:`FonPaqLop06:hypervolume` had the HV3D algorithm as the base + case, giving a time complexity of :math:`O(n^{m-2} \log n)`. Andreia + P. Guerreiro enhanced the numerical stability of the recursive algorithm by + avoiding floating-point comparisons of partial hypervolumes. References diff --git a/r/NEWS.md b/r/NEWS.md index c5c5e8868..7146c15ef 100644 --- a/r/NEWS.md +++ b/r/NEWS.md @@ -7,6 +7,7 @@ * Fix documentation of `epsilon_additive()` (@leandrolanzieri) * `hypervolume()` uses the inclusion-exclusion algorithm for small inputs of up to 15 points, which is significantly faster. + * `hypervolume()` is significantly faster for more than four dimensions (Andreia P. Guerreiro). # moocore 0.1.10 diff --git a/r/R/hv.R b/r/R/hv.R index 85651be23..9cfd69874 100644 --- a/r/R/hv.R +++ b/r/R/hv.R @@ -61,15 +61,15 @@ #' inclusion-exclusion algorithm \citep{WuAzam2001metrics}, which has \eqn{O(m #' 2^{n})} time and \eqn{O(n\cdot m)} space complexity, but it is very fast for #' such small sets. For larger number of points, it uses a recursive algorithm -#' \citep{FonPaqLop06:hypervolume} with \eqn{\text{HV4D}^{+}} as the base case, -#' resulting in a \eqn{O(n^{m-2})} time complexity and \eqn{O(n)} space -#' complexity in the worst-case. Experimental results show that the pruning -#' techniques used may reduce the time complexity even further. The original -#' proposal \citep{FonPaqLop06:hypervolume} had the HV3D algorithm as the base -#' case, giving a time complexity of \eqn{O(n^{m-2} \log n)}. Andreia -#' P. Guerreiro enhanced the numerical stability of the algorithm by avoiding -#' floating-point comparisons of partial hypervolumes. -#' +#' \citep{FonPaqLop06:hypervolume} that computes 4D contributions +#' \citep{GueFon2017hv4d} as the base case, resulting in a \eqn{O(n^{m-2})} +#' time complexity and \eqn{O(n)} space complexity in the worst-case. +#' Experimental results show that the pruning techniques used may reduce the +#' time complexity even further. The original proposal +#' \citep{FonPaqLop06:hypervolume} had the HV3D algorithm as the base case, +#' giving a time complexity of \eqn{O(n^{m-2} \log n)}. Andreia P. Guerreiro +#' enhanced the numerical stability of the algorithm by avoiding floating-point +#' comparisons of partial hypervolumes. #' #' @references #' diff --git a/r/man/hypervolume.Rd b/r/man/hypervolume.Rd index 36e8f8a18..8113b6f02 100644 --- a/r/man/hypervolume.Rd +++ b/r/man/hypervolume.Rd @@ -72,14 +72,15 @@ For 5D or higher and up to 15 points, the implementation uses the inclusion-exclusion algorithm \citep{WuAzam2001metrics}, which has \eqn{O(m 2^{n})} time and \eqn{O(n\cdot m)} space complexity, but it is very fast for such small sets. For larger number of points, it uses a recursive algorithm -\citep{FonPaqLop06:hypervolume} with \eqn{\text{HV4D}^{+}} as the base case, -resulting in a \eqn{O(n^{m-2})} time complexity and \eqn{O(n)} space -complexity in the worst-case. Experimental results show that the pruning -techniques used may reduce the time complexity even further. The original -proposal \citep{FonPaqLop06:hypervolume} had the HV3D algorithm as the base -case, giving a time complexity of \eqn{O(n^{m-2} \log n)}. Andreia -P. Guerreiro enhanced the numerical stability of the algorithm by avoiding -floating-point comparisons of partial hypervolumes. +\citep{FonPaqLop06:hypervolume} that computes 4D contributions +\citep{GueFon2017hv4d} as the base case, resulting in a \eqn{O(n^{m-2})} +time complexity and \eqn{O(n)} space complexity in the worst-case. +Experimental results show that the pruning techniques used may reduce the +time complexity even further. The original proposal +\citep{FonPaqLop06:hypervolume} had the HV3D algorithm as the base case, +giving a time complexity of \eqn{O(n^{m-2} \log n)}. Andreia P. Guerreiro +enhanced the numerical stability of the algorithm by avoiding floating-point +comparisons of partial hypervolumes. } \examples{