From 9b33f37aa44e40d3dab0a820a64ac5bc7403c41a Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Tue, 25 Nov 2025 17:02:02 +0000 Subject: [PATCH 1/6] Compute a 4D hv contribution in the base case. Improve the 4D base case of FPL: * c/hv_priv.h: Remove x_aux from dlnote_t, x is const again. * c/hv.c: - Remove references to x_aux. - Also maintain a list sorted according to z coordinate. - Allocate (once) the auxiliar data structures used in the base case. * c/hv4d_priv.h (onec4dplusU): - Reconstruct the 3D base of the 4D contribution by sweeping the (sorted) z list. - Use an auxiliar dlnode_t list (list_aux) and an auxiliar double vector (x_aux) for modified points. - Take advantage of sweeping points in ascending order of the z-coordinate (continue_base_update_z_closest). * c/hv.c: Use ->vol instead of ->x for the auxiliary vectors. (fpli_hv4d): Delete unused. * c/hv4d_priv.h (update_bound_3d): New. --- c/hv.c | 130 +++++++++----- c/hv4d_priv.h | 316 ++++++++++++++++++++++++++++++++- c/sort.h | 6 + python/src/moocore/_moocore.py | 2 +- r/R/hv.R | 14 +- r/man/hypervolume.Rd | 13 +- 6 files changed, 415 insertions(+), 66 deletions(-) diff --git a/c/hv.c b/c/hv.c index 2c3fc64f..248962f1 100644 --- a/c/hv.c +++ b/c/hv.c @@ -80,7 +80,6 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, if (unlikely(n <= MAX_ROWS_HV_INEX)) goto finish; - // Allocate single blocks of memory as much as possible. // 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)); @@ -96,7 +95,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. + + // Reserve space for auxiliar list and auxiliar x vector used in onec4dplusU. + // This list will store at most n-1 points. + dlnode_t * list_aux = (dlnode_t *) malloc(n * sizeof(*list_aux)); + list_aux->vol = malloc(3 * n * sizeof(double)); + 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. @@ -120,25 +125,28 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, */ // 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; + if (j <= -1) { + const int j2 = j + 2; + (list4d+1)->next[j2] = scratch[0]; + scratch[0]->prev[j2] = list4d+1; for (i = 1; i < n; i++) { - scratch[i-1]->next[1] = scratch[i]; - scratch[i]->prev[1] = scratch[i-1]; + scratch[i-1]->next[j2] = scratch[i]; + scratch[i]->prev[j2] = 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++) { - scratch[i-1]->r_next[j] = scratch[i]; - scratch[i]->r_prev[j] = scratch[i-1]; + scratch[n-1]->next[j2] = list4d+2; + (list4d+2)->prev[j2] = scratch[n-1]; + if (j == -2) + break; + } else { + head->r_next[j] = scratch[0]; + scratch[0]->r_prev[j] = head; + for (i = 1; i < n; i++) { + scratch[i-1]->r_next[j] = scratch[i]; + scratch[i]->r_prev[j] = scratch[i-1]; + } + scratch[n-1]->r_next[j] = head; + head->r_prev[j] = scratch[n-1]; } - 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++) @@ -146,7 +154,7 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } // Reset x to point to the first objective. for (i = 1; i <= n; i++) - head[i].x -= STOP_DIMENSION; + head[i].x -= STOP_DIMENSION - 1; free(scratch); @@ -161,8 +169,10 @@ 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->prev[0]->vol); // free x_aux (4D basecase) + free(head->prev[0]); // free list_aux (4D basecase) free(head->r_next); free(head->area); free(head); @@ -180,6 +190,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 +228,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 +250,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 +262,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 @@ -357,7 +392,7 @@ _attr_maybe_unused static size_t debug_counter[6] = { 0 }; static double hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, - const double * restrict ref, double * restrict bound) + const double * restrict ref, double * restrict bound, dlnode_t * the_point) { ASSUME(c > 1); ASSUME(dim >= STOP_DIMENSION); @@ -365,7 +400,7 @@ hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, /*--------------------------------------- base case of dimension 4 --------------------------------------*/ - return fpli_hv4d(list, c); + return fpli_onec4d(list, c, the_point); } ASSUME(dim > STOP_DIMENSION); /* ------------------------------------------------------ @@ -434,7 +469,10 @@ 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); + hypera = hv_recursive(list, dim - 1, c, ref, bound, p1); + if (dim - 1 == STOP_DIMENSION) { // hypera only has the contribution of p1 + hypera += p1_prev->area[d_stop]; + } /* 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. */ @@ -495,19 +533,19 @@ fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, // FIXME: This is never used? // bound[d_stop] = p0->x[dim]; reinsert_nobound(p0, dim); - p1 = p0; - dlnode_t * p1_prev = p0->r_prev[d_stop - 1]; - p0 = p0->r_next[d_stop - 1]; - c++; - assert(c > 1); while (true) { - // FIXME: This is not true in the first iteration if c > 1 previously. - //assert(p0 == p1->r_prev[d_stop]); - assert(p1_prev == p1->r_prev[d_stop - 1]); + dlnode_t * p1_prev = p0->r_prev[d_stop - 1]; + p1 = p0; + p0 = p0->r_next[d_stop - 1]; p1->vol[d_stop] = hyperv; assert(p1->ignore == 0); - double hypera = hv_recursive(list, dim - 1, c, ref, bound); + c++; + + double hypera = hv_recursive(list, dim - 1, c, ref, bound, p1); + if (dim - 1 == STOP_DIMENSION) { //hypera only has the contribution of p1 + hypera += p1_prev->area[d_stop]; + } /* 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. */ @@ -533,10 +571,6 @@ fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, // bound[d_stop] = p0->x[dim]; // FIXME: Does updating the bound here matters? reinsert(p0, dim, bound); - p1 = p0; - p1_prev = p0->r_prev[d_stop - 1]; - p0 = p0->r_next[d_stop - 1]; - c++; } } diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 0adaddfc..1bf381d7 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -5,7 +5,7 @@ HV4D+ algorithm. ------------------------------------------------------------------------------ - Copyright (C) 2013, 2016, 2017 + Copyright (C) 2013, 2016, 2017, 2025 Andreia P. Guerreiro This Source Code Form is subject to the terms of the Mozilla Public @@ -69,7 +69,7 @@ update_links(dlnode_t * restrict list, dlnode_t * restrict newp) assert(list+2 == list->prev[0]); const double newx[] = { newp->x[0], newp->x[1], newp->x[2] }; dlnode_t * p = newp->next[0]; - const dlnode_t * stop = list+2; + const dlnode_t * const stop = list+2; while (p != stop) { const double * px = p->x; // px dominates newx (but not equal) @@ -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,12 +130,18 @@ 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; } //if (!lexicographic_less_3d(px, newx)) { // Slower if (!(p_lt_new_2 || (p_eq_new_2 && (p_lt_new_1 || (p_eq_new_1 && p_leq_new_0))))) { + assert(!lexicographic_less_3d(px, newx)); newp->closest[0] = closest0; newp->closest[1] = closest1; newp->prev[0] = p->prev[0]; @@ -159,6 +169,91 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n unreachable(); } + +/** + 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; +} + + // FIXME: This is very similar to the loop in hvc3d_list() but it doesn't use p->last_slice_z static double one_contribution_3d(dlnode_t * restrict newp) @@ -216,7 +311,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]); @@ -225,7 +320,7 @@ hv4dplusU(dlnode_t * list) double volume = 0, hv = 0; dlnode_t * newp = (list+1)->next[1]; - const dlnode_t * last = list+2; + const dlnode_t * const last = list+2; while (newp != last) { if (restart_base_setup_z_and_closest(list, newp)) { // newp was not dominated by something else. @@ -245,4 +340,217 @@ hv4dplusU(dlnode_t * list) return hv; } +static inline int compare_lex_2d(const void * pa, const void * pb) +{ + const double ax = **(const double **)pa; + const double bx = **(const double **)pb; + const double ay = *(*(const double **)pa + 1); + const double by = *(*(const double **)pb + 1); + int cmpx = cmp_double_asc(ax, bx); + int cmpy = cmp_double_asc(ay, by); + return cmpy ? cmpy : cmpx; +} + +static inline void lex_sort_equal_z_and_setup_nodes(dlnode_t * newp_aux, double * x_aux, size_t n) { + const double **scratch = malloc(n * sizeof(*scratch)); + double * x = x_aux; + size_t i; + for (i = 0; i < n; i++){ + scratch[i] = x; + x += 3; + } + + qsort(scratch, n, sizeof(*scratch), compare_lex_2d); + + for (i = 0; i < n; i++) { + newp_aux->x = scratch[i]; + assert(i == 0 || lexicographic_less_2d(scratch[i-1], scratch[i])); + newp_aux++; + } + + free(scratch); +} + + +static inline void +update_bound_3d(double * restrict dest, const double * restrict a, const double * restrict b) +{ + for (int i = 0; i < 3; i++) + dest[i] = MAX(a[i], b[i]); +} + +// FIXME: Move this function to hv.c +#ifdef HV_RECURSIVE +/** + 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]); + + dlnode_t * newp = (list+1)->next[0]; + const dlnode_t * const last = list+2; + dlnode_t * const z_first = newp; + 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 + dlnode_t * prevp_aux; + int c; + + 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) { + 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; + 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; + } + + // FIXME: Is it possible that x_aux matches newpx? YES! So just assign and avoid the comparison. + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + update_bound_3d(x_aux, newpx, the_point_x); + 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 + c = 0; + while (newp != last) { + const double * newpx = newp->x; + if (newpx[3] <= the_point_x[3] && newp->ignore < 3) { + + // FIXME: Is it possible that x_aux matches newpx? YES! So just assign and avoid the comparison. + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + update_bound_3d(x_aux, newpx, the_point_x); + 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(newp_aux, x_aux-3*c, c); + continue_base_update_z_closest(list, newp_aux, false); + prevp_aux = newp_aux; + 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 = newp_aux; + newp_aux++; + } + } + c = 0; + } + newp = newp->next[0]; + } + } + + 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]; + // FIXME: Does this call always return true? +#if DEBUG >= 1 + assert(restart_base_setup_z_and_closest(list, the_point)); +#else + restart_base_setup_z_and_closest(list, the_point); +#endif + 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) { + // FIXME: Is it possible that x_aux matches newpx? YES! So just assign and avoid the comparison. + // x_aux is the coordinate-wise maximum between newpx and the_point_x. + update_bound_3d(x_aux, newpx, the_point_x); + 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++; + } + } + // 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 // HV_RECURSIVE + #endif // _HV4D_PRIV_H diff --git a/c/sort.h b/c/sort.h index 6c9221fb..32631f1e 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) { diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index 1c672b03..a7c76405 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -465,7 +465,7 @@ 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, + :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 diff --git a/r/R/hv.R b/r/R/hv.R index 85651be2..1c3698ae 100644 --- a/r/R/hv.R +++ b/r/R/hv.R @@ -61,16 +61,16 @@ #' 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 +#' \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 #' #' \insertAllCited{} diff --git a/r/man/hypervolume.Rd b/r/man/hypervolume.Rd index 36e8f8a1..c6129476 100644 --- a/r/man/hypervolume.Rd +++ b/r/man/hypervolume.Rd @@ -72,12 +72,13 @@ 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 +\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. } From 7669e76286cc0732b49b09bb6b23f561b55b0ad2 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Thu, 18 Dec 2025 00:10:38 +0000 Subject: [PATCH 2/6] hv4d_priv.h: Detect dominated points and set ignore to 3. --- c/hv4d_priv.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 1bf381d7..94b098a8 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -371,7 +371,6 @@ static inline void lex_sort_equal_z_and_setup_nodes(dlnode_t * newp_aux, double free(scratch); } - static inline void update_bound_3d(double * restrict dest, const double * restrict a, const double * restrict b) { @@ -535,7 +534,12 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_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]; From 027d3e532ad62e92abfe4123cdc634088955b4bb Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 26 Dec 2025 10:45:49 +0000 Subject: [PATCH 3/6] debugging: use -O0 --- python/src/moocore/_ffi_build.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/src/moocore/_ffi_build.py b/python/src/moocore/_ffi_build.py index 9c6ea66a..4a560cc4 100644 --- a/python/src/moocore/_ffi_build.py +++ b/python/src/moocore/_ffi_build.py @@ -85,7 +85,7 @@ def _get_target_platform(): # Compiler flags. MSVC_CFLAGS = ["/GL", "/O2", "/GS-", "/wd4996", "/DMOOCORE_SHARED_LIB"] MSVC_LDFLAGS = ["/LTCG"] -GCC_CFLAGS = ["-O3", "-flto", "-fvisibility=hidden"] +GCC_CFLAGS = ["-O0", "-flto", "-fvisibility=hidden"] GCC_LDFLAGS = [*GCC_CFLAGS] if is_x86_64: # Compile for sufficiently old x86-64 architecture. From 5feb5e373c9829e26451b72f1fe9af2318f9045c Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 26 Dec 2025 11:10:41 +0000 Subject: [PATCH 4/6] Disable hv-inex for debugging. --- c/hv.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/c/hv.c b/c/hv.c index 248962f1..0c7c4764 100644 --- a/c/hv.c +++ b/c/hv.c @@ -37,8 +37,8 @@ #define HV_RECURSIVE #include "hv4d_priv.h" -#define STOP_DIMENSION 3 // default: stop on dimension 4. -#define MAX_ROWS_HV_INEX 15 +#define STOP_DIMENSION 3 // stop on dimension 4. +#define MAX_ROWS_HV_INEX 1 static int compare_node(const void * restrict p1, const void * restrict p2) { From 0819cf8536cfd3c3d626457a73c2b37a108df727 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 26 Dec 2025 11:38:14 +0000 Subject: [PATCH 5/6] enable asan debugging --- c/Makefile | 22 ++++++++++++++++++++++ c/gcc.mk | 2 +- c/hv.c | 4 ++-- 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/c/Makefile b/c/Makefile index f31e54a6..87ec76fa 100644 --- a/c/Makefile +++ b/c/Makefile @@ -156,6 +156,7 @@ EXE_CFLAGS += $(SANITIZERS) $(OPT_CFLAGS) $(MARCH_FLAGS) $(WARN_CFLAGS) \ # List of test target base names TEST_NAMES := eaf nondominated hv hvapprox igd epsilon dominatedsets ndsort +ASAN_TARGETS := $(addprefix asan-,$(TEST_NAMES)) TEST_TARGETS := $(addprefix test-,$(TEST_NAMES)) TIME_TARGETS := $(addprefix time-,$(TEST_NAMES)) ALL_NAMES := $(TEST_NAMES) @@ -242,6 +243,27 @@ clean: @$(RM) $(OBJS) $(HV_OBJS) $(HV_LIB) TESTSUITE:="../testsuite" +asan: OPT_CFLAGS ?=-Og -g3 +asan: DEBUG=1 +asan: SANITIZERS ?=-fsanitize=undefined,address,pointer-compare,pointer-subtract,float-cast-overflow,float-divide-by-zero +asan: + @if ! test -d $(TESTSUITE); then \ + echo "Error: Testsuite not found in $(TESTSUITE)" && exit 1; \ + fi + $(MAKE) all DEBUG=1 SANITIZERS="$(SANITIZERS)" OPT_CFLAGS="$(OPT_CFLAGS)" + $(MAKE) $(ASAN_TARGETS) + +$(ASAN_TARGETS): target=$(subst asan-,,$@) +$(ASAN_TARGETS): DEBUG=1 +$(ASAN_TARGETS): OPT_CFLAGS ?=-Og -g3 +$(ASAN_TARGETS): SANITIZERS ?=-fsanitize=undefined,address,pointer-compare,pointer-subtract,float-cast-overflow,float-divide-by-zero +$(ASAN_TARGETS): + $(MAKE) $(target) DEBUG=1 SANITIZERS="$(SANITIZERS)" OPT_CFLAGS="$(OPT_CFLAGS)" + $(call ECHO,---> Testing $(target) for errors (DEBUG=1, $(SANITIZERS)) <--- ) + export ASAN_OPTIONS='detect_invalid_pointer_pairs=2:abort_on_error=1' + export UBSAN_OPTIONS='print_stacktrace=1:halt_on_error=1' + @cd $(TESTSUITE)/$(target) && ../regtest.py $(BINDIR)/$(target)$(EXE) && cd $(PWD) + test: OPT_CFLAGS ?=-Og -g3 test: SANITIZERS ?= test: diff --git a/c/gcc.mk b/c/gcc.mk index e68f3da0..357a056b 100644 --- a/c/gcc.mk +++ b/c/gcc.mk @@ -13,7 +13,7 @@ ifeq ($(DEBUG), 0) # -ftree-loop-im -ftree-loop-ivcanon -fivopts -ftree-vectorize -funroll-loops -fipa-pta # else - SANITIZERS ?= -fsanitize=undefined -fsanitize=address -fsanitize=float-cast-overflow -fsanitize=float-divide-by-zero + SANITIZERS ?= -fsanitize=undefined,address,pointer-compare,pointer-subtract,float-cast-overflow,float-divide-by-zero OPT_CFLAGS ?= -g3 -O0 endif diff --git a/c/hv.c b/c/hv.c index 0c7c4764..7c5e40ae 100644 --- a/c/hv.c +++ b/c/hv.c @@ -38,7 +38,7 @@ #include "hv4d_priv.h" #define STOP_DIMENSION 3 // stop on dimension 4. -#define MAX_ROWS_HV_INEX 1 +#define MAX_ROWS_HV_INEX 2 static int compare_node(const void * restrict p1, const void * restrict p2) { @@ -356,7 +356,7 @@ hv_inex_list(const dlnode_t * restrict list, int n, dimension_t dim, child = subset_max[top]; upper_bound(child, list[idx].x, parent, dim); // Inclusion–exclusion accumulation. - hv[top & 1] += one_point_hv(child, ref, dim); + hv[top % 2] += one_point_hv(child, ref, dim); idx++; } else if (top > 0) { --top; From 74063d75b791c06aa907efd4a02df903cd7c3d2d Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 26 Dec 2025 16:00:10 +0000 Subject: [PATCH 6/6] debugging: remove nondominated --- c/hv.c | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/c/hv.c b/c/hv.c index 7c5e40ae..d33c16d8 100644 --- a/c/hv.c +++ b/c/hv.c @@ -34,6 +34,7 @@ #include #include "common.h" #include "hv.h" +#include "nondominated.h" #define HV_RECURSIVE #include "hv4d_priv.h" @@ -64,6 +65,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, dimension_t d_stop = d - STOP_DIMENSION; size_t n = *size; + bool * nondom = is_nondominated_minimise(data, n, d, /* keep_weakly=*/false); + dlnode_t * head = malloc((n+1) * sizeof(*head)); size_t i = 1; for (size_t j = 0; j < n; j++) { @@ -71,11 +74,12 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, point. This is needed to assure that the points left are only those that are needed to calculate the hypervolume. */ const double * restrict px = data + j * d; - if (likely(strongly_dominates(px, ref, d))) { + if (nondom[j] && likely(strongly_dominates(px, ref, d))) { head[i].x = px; // this will be fixed a few lines below... i++; } } + free(nondom); n = i - 1; if (unlikely(n <= MAX_ROWS_HV_INEX)) goto finish;