From 28a2bed648fab645c292acdd7fdf802c4c4b4794 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 28 Dec 2025 09:03:10 +0000 Subject: [PATCH 01/19] * c/hv.c (hv_recursive): Do not recurse directly into basecase. --- c/hv.c | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/c/hv.c b/c/hv.c index ac043e285..8b7186451 100644 --- a/c/hv.c +++ b/c/hv.c @@ -360,14 +360,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 +426,13 @@ 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_hv4d(list, c); + } 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 +501,14 @@ 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_hv4d(list, c); + } 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. */ From a7376cad6b186b2cc93ac80854e1529fa1b94b86 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Tue, 25 Nov 2025 17:02:02 +0000 Subject: [PATCH 02/19] 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: Simplify code, fix warnings, add more asserts and comments. (update_bound_3d): New. --- c/hv.c | 111 ++++++++++++++++-------- c/hv4d_priv.h | 231 +++++++++++++++++++++++++++++++++++++++++++++++++- c/sort.h | 6 ++ 3 files changed, 310 insertions(+), 38 deletions(-) diff --git a/c/hv.c b/c/hv.c index 8b7186451..e9527b669 100644 --- a/c/hv.c +++ b/c/hv.c @@ -96,7 +96,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. @@ -112,41 +118,41 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, 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 >= -2; 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; + if (j >= 0) { + head->r_next[j] = scratch[0]; + scratch[0]->r_prev[j] = head; for (i = 1; i < n; i++) { - scratch[i-1]->next[1] = scratch[i]; - scratch[i]->prev[1] = scratch[i-1]; + scratch[i-1]->r_next[j] = scratch[i]; + scratch[i]->r_prev[j] = 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]->r_next[j] = head; + head->r_prev[j] = scratch[n-1]; + } else { + 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[j2] = scratch[i]; + scratch[i]->prev[j2] = scratch[i-1]; + } + scratch[n-1]->next[j2] = list4d+2; + (list4d+2)->prev[j2] = 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++) head[i].x--; } // Reset x to point to the first objective. - for (i = 1; i <= n; i++) - head[i].x -= STOP_DIMENSION; + for (i = 1; i <= n; i++) { + head[i].x -= STOP_DIMENSION - 2; + } free(scratch); @@ -161,8 +167,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 +188,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 +226,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 +248,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 +260,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 @@ -429,7 +462,9 @@ hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, ASSUME(dim - 1 >= STOP_DIMENSION); if (dim - 1 == STOP_DIMENSION) { // base case of dimension 4. - hypera = fpli_hv4d(list, c); + 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); } @@ -505,7 +540,9 @@ fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, ASSUME(dim - 1 >= STOP_DIMENSION); if (dim - 1 == STOP_DIMENSION) { // base case of dimension 4. - hypera = fpli_hv4d(list, c); + 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); } diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index f92b8a340..fc58536ba 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; } @@ -162,6 +171,96 @@ 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 * lex_prev = list+1; + dlnode_t * p = lex_prev->cnext[1]; + dlnode_t * closest1; + + // 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)); + + lex_prev = p; // newp->closest[0] = lex_prev + + // Check if newp is dominated. + if (weakly_dominates(p->x, newx, 3)) { + return false; + } + + 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]; + + closest1 = 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]); + p = p->cnext[1]; + } + closest1 = p; + } + + newp->closest[0] = lex_prev; + newp->closest[1] = closest1; + + // 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) @@ -219,7 +318,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]); @@ -248,4 +347,134 @@ hv4dplusU(dlnode_t * list) return hv; } + +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]; + dlnode_t * z_first = newp; + const dlnode_t * const last = list+2; + dlnode_t * 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); + + // FIXME: Would it be possible to remove the_point so we don't need to test it? + + const double * the_point_x = the_point->x; + if ((list+1)->next[1] != the_point) { + bool done_once = false; + // PART 1: Setup 3D base (if there are any points below the_point_x[3]) + assert(newp != last); + do { + const double * newpx = newp->x; + if (newpx[3] <= the_point_x[3] && newp != the_point && newp->ignore < 3){ + if (newpx[0] <= the_point_x[0] && newpx[1] <= the_point_x[1] && newpx[2] <= the_point_x[2]) { + 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. + update_bound_3d(x_aux, newpx, the_point_x); + newp_aux->x = x_aux; + x_aux += 3; + bool equalz = newp_aux->x[2] == the_point_x[2] && done_once; + if (continue_base_update_z_closest(list, newp_aux, equalz)) { + newp_aux++; + done_once = true; + } + } + newp = newp->next[0]; + } while (newp != last); + } + 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: Why do we ignore the return value of this function? + bool res = restart_base_setup_z_and_closest(list, the_point); + assert(res); + 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 2: Update the 3D contribution. + while (newp != last) { + const double * newpx = newp->x; + if (newpx[0] <= the_point_x[0] && newpx[1] <= the_point_x[1] && newpx[2] <= the_point_x[2]) + break; + + if (newp->ignore < 3) { + // 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); + } + // FIXME: If the above returned false, then newp_aux was not used + // so we could re-use it, no? + 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 3f78fd5a7..ca211d5b4 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) { From fcf8154856cde5f6053f2f31319ebfb889a71c73 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 13 Dec 2025 13:09:21 +0000 Subject: [PATCH 03/19] Add documentation. --- python/src/moocore/_moocore.py | 17 +++++++++-------- r/R/hv.R | 18 +++++++++--------- r/man/hypervolume.Rd | 17 +++++++++-------- 3 files changed, 27 insertions(+), 25 deletions(-) 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/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{ From e0be01af6e33d8d1bc0a87032b0b9410d6854f62 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 13 Dec 2025 16:24:30 +0000 Subject: [PATCH 04/19] more simplifications avoid newp != the_point. more cleanups hv4d_priv.h: Move newp_aux++ within the if(). hv4d_priv.h: Avoid warning with DEBUG=0. --- c/hv4d_priv.h | 62 +++++++++++++++++++++++++-------------------------- 1 file changed, 30 insertions(+), 32 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index fc58536ba..71ef131c6 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -182,11 +182,9 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n 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]}; + const double newx[] = { newp->x[0], newp->x[1], newp->x[2] }; assert(list+1 == list->next[0]); - dlnode_t * lex_prev = list+1; - dlnode_t * p = lex_prev->cnext[1]; - dlnode_t * closest1; + dlnode_t * p = (list+1)->cnext[1]; // find newp->closest[0] while (p->x[1] < newx[1]){ @@ -200,13 +198,11 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new assert(lexicographic_less_2d(p->x, newx)); assert(lexicographic_less_2d(newx, p->cnext[1]->x)); - lex_prev = p; // newp->closest[0] = lex_prev - // 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] @@ -229,24 +225,22 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new newp->prev[0] = lex_prev; newp->next[0] = lex_prev->next[0]; - closest1 = list; + 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]); - p = p->cnext[1]; + // find newp->closest[1] + while (p->x[0] >= newx[0]){ + assert(p->x[2] <= newx[2]); + p = p->cnext[1]; } - closest1 = p; + newp->closest[1] = p; } newp->closest[0] = lex_prev; - newp->closest[1] = closest1; - // update cnext lex_prev->cnext[1] = newp; newp->cnext[0] = lex_prev; @@ -375,10 +369,9 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, assert(list+1 == list->next[1]); dlnode_t * newp = (list+1)->next[0]; - dlnode_t * z_first = newp; const dlnode_t * const last = list+2; - dlnode_t * z_last = last->prev[0]; - + 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 @@ -386,24 +379,26 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, reset_sentinels_3d(list); restart_list_y(list); - // FIXME: Would it be possible to remove the_point so we don't need to test it? - const double * the_point_x = the_point->x; if ((list+1)->next[1] != the_point) { - bool done_once = false; // PART 1: Setup 3D base (if there are any points below the_point_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; assert(newp != last); do { const double * newpx = newp->x; - if (newpx[3] <= the_point_x[3] && newp != the_point && newp->ignore < 3){ - if (newpx[0] <= the_point_x[0] && newpx[1] <= the_point_x[1] && newpx[2] <= the_point_x[2]) { - the_point->ignore = 3; + 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; @@ -416,6 +411,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, } newp = newp->next[0]; } while (newp != last); + the_point->ignore = the_point_ignore; } newp = the_point->next[1]; while (newp->x[3] <= the_point_x[3]) @@ -423,9 +419,12 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, dlnode_t * tp_prev_z = the_point->prev[0]; dlnode_t * tp_next_z = the_point->next[0]; - // FIXME: Why do we ignore the return value of this function? - bool res = restart_base_setup_z_and_closest(list, the_point); - assert(res); + // 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; @@ -439,10 +438,11 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, // PART 2: Update the 3D contribution. while (newp != last) { const double * newpx = newp->x; - if (newpx[0] <= the_point_x[0] && newpx[1] <= the_point_x[1] && newpx[2] <= the_point_x[2]) + 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; @@ -455,10 +455,8 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, add_to_z(newp_aux); update_links(list, newp_aux); + newp_aux++; } - // FIXME: If the above returned false, then newp_aux was not used - // so we could re-use it, no? - newp_aux++; } // FIXME: If newp was dominated, can we accumulate the height and update // hv later? AG: Yes From 26570cfdb5a4c965cb7f96ea44c342794a444936 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Wed, 17 Dec 2025 16:06:47 +0000 Subject: [PATCH 05/19] * c/hv4d_priv.h: Fix bug related to repeated coordinates. --- c/hv4d_priv.h | 92 +++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 86 insertions(+), 6 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 71ef131c6..94b952eae 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -235,6 +235,7 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new // 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; @@ -341,6 +342,37 @@ 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) @@ -375,20 +407,25 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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) { - // PART 1: Setup 3D base (if there are any points below the_point_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; assert(newp != last); - do { + + // 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)) { @@ -398,21 +435,64 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, (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; - bool equalz = newp_aux->x[2] == the_point_x[2] && done_once; - if (continue_base_update_z_closest(list, newp_aux, equalz)) { + + if (continue_base_update_z_closest(list, newp_aux, done_once)) { newp_aux++; done_once = true; } } newp = newp->next[0]; - } while (newp != last); + } 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]; @@ -435,7 +515,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, assert(height > 0); double hv = volume * height; - // PART 2: Update the 3D contribution. + // PART 3: Update the 3D contribution. while (newp != last) { const double * newpx = newp->x; if (weakly_dominates(newpx, the_point_x, 3)) From a02658fd626daea3590a12cc8f13126f65bba1f2 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Thu, 18 Dec 2025 00:10:38 +0000 Subject: [PATCH 06/19] 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 94b952eae..6ab231bec 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -373,7 +373,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) { @@ -537,7 +536,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 fd8323827f45202760f937d4a88f70788627f01b Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 29 Dec 2025 19:48:03 +0000 Subject: [PATCH 07/19] * c/hv.c (fpli_setup_cdllist): Split for-loop for clarity. --- c/hv.c | 53 +++++++++++++++++++++++++++++------------------------ 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/c/hv.c b/c/hv.c index e9527b669..0ea48aa8a 100644 --- a/c/hv.c +++ b/c/hv.c @@ -118,42 +118,47 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, for (i = 0; i < n; i++) scratch[i] = head + 1 + i; - for (int j = d_stop - 2; j >= -2; j--) { + 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 >= 0) { - 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]; - } else { - 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[j2] = scratch[i]; - scratch[i]->prev[j2] = scratch[i-1]; - } - scratch[n-1]->next[j2] = list4d+2; - (list4d+2)->prev[j2] = scratch[n-1]; + 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]; // 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 - 2; - } + 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. From b00c1d08528ac62a52e9a7911ca1fe964b0ba953 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 7 Jan 2026 15:07:57 +0000 Subject: [PATCH 08/19] bugfix? --- c/hv4d_priv.h | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 6ab231bec..b8c27e51c 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -267,7 +267,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); @@ -399,9 +403,8 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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_first = (list+1)->next[0]; dlnode_t * const z_last = last->prev[0]; double * x_aux = list_aux->vol; @@ -412,15 +415,15 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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) { + if ((list+1)->next[1]->x[3] <= the_point_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 @@ -434,8 +437,6 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, (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; @@ -455,17 +456,15 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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 > 0 && newp->next[0]->x[2] > newpx[2]) { if (c == 1) { - newp_aux->x = x_aux-3; + newp_aux->x = x_aux - 3; continue_base_update_z_closest(list, newp_aux, false); newp_aux++; } else { @@ -492,7 +491,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, } } - newp = the_point->next[1]; + dlnode_t * newp = the_point->next[1]; while (newp->x[3] <= the_point_x[3]) newp = newp->next[1]; From d9096f2d8272a74c32f0e8695239b0f1554dc68a Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Wed, 7 Jan 2026 15:09:04 +0000 Subject: [PATCH 09/19] *c/hv4d_priv.h: Fix bug related to repeated coordinates --- c/hv4d_priv.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index b8c27e51c..56f875d4c 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -417,7 +417,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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]->x[3] <= 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. @@ -518,7 +518,6 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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. From 0be2fe1650096eb97f056806877ececd42ae127e Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 7 Jan 2026 17:31:27 +0000 Subject: [PATCH 10/19] more minor cleanups --- c/hv4d_priv.h | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 56f875d4c..cf4784694 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -409,14 +409,12 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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]) + // 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 @@ -452,7 +450,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, the_point->ignore = the_point_ignore; // PART 2: Setup the remainder of the 3D base - c = 0; + int c = 0; while (newp != last) { const double * newpx = newp->x; if (newpx[3] <= the_point_x[3] && newp->ignore < 3) { @@ -469,19 +467,19 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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); + 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; + 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]){ + 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; + prevp_aux_x = newp_aux->x; newp_aux++; } } @@ -497,7 +495,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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? + // This call should always return true. #if DEBUG >= 1 assert(restart_base_setup_z_and_closest(list, the_point)); #else @@ -519,7 +517,6 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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; From 570cf8bd858e99f4f8e7507888ae7522d1940c6e Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 15:12:12 +0000 Subject: [PATCH 11/19] c/NEWS.md: Document faster fpli_hv(). --- c/NEWS.md | 2 ++ 1 file changed, 2 insertions(+) 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 From 21e5744a971b50703a2fba1c1df0d75c87f93d69 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 15:37:42 +0000 Subject: [PATCH 12/19] c/hv.c (fpli_setup_cdllist): Do not allocate memory for an extra point. --- c/hv.c | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/c/hv.c b/c/hv.c index 0ea48aa8a..b486b47bd 100644 --- a/c/hv.c +++ b/c/hv.c @@ -86,8 +86,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, 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); + head->area = malloc(2 * d_stop * n * sizeof(*data)); + head->vol = head->area + d_stop * n; head->x = NULL; // head contains no data head->ignore = 0; // should never get used @@ -110,9 +110,12 @@ 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++) @@ -161,10 +164,6 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } 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; @@ -177,7 +176,7 @@ static void fpli_free_cdllist(dlnode_t * head) 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[1].area); // free ->area and ->vol. free(head); } From 3a97a9e0f59da8281553c1e2e8f4c19fe21045e0 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 15:49:45 +0000 Subject: [PATCH 13/19] c/hv.c (fpli_setup_cdllist): Allocate memory for main and auxiliary CDLLs together. --- c/hv.c | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/c/hv.c b/c/hv.c index b486b47bd..b4fe8a2cb 100644 --- a/c/hv.c +++ b/c/hv.c @@ -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 @@ -97,9 +101,9 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, // should remain untouched. head->next[0] = list4d; - // 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)); + // Setup the auxiliary list used in onec4dplusU(). + dlnode_t * list_aux = head + n + 1; + // Reserve space for n auxiliary 3D points used by onec4dplusU(). list_aux->vol = malloc(3 * n * sizeof(double)); head->prev[0] = list_aux; list_aux->next[0] = list4d; @@ -174,10 +178,9 @@ static void fpli_free_cdllist(dlnode_t * head) 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[1].area); // free ->area and ->vol. - free(head); + free(head); // Free main CDLL and list_aux (4D basecase). } static inline void From d45d9e2dc7180cfc20b2286a3cca93b6d463fc56 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 16:06:42 +0000 Subject: [PATCH 14/19] c/hv.c (fpli_setup_cdllist): Allocate together the memory for ->area, ->vol and auxiliary points. --- c/hv.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/c/hv.c b/c/hv.c index b4fe8a2cb..2f36f0f80 100644 --- a/c/hv.c +++ b/c/hv.c @@ -89,8 +89,9 @@ 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 * sizeof(*data)); + // 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 @@ -103,8 +104,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, // Setup the auxiliary list used in onec4dplusU(). dlnode_t * list_aux = head + n + 1; - // Reserve space for n auxiliary 3D points used by onec4dplusU(). - list_aux->vol = malloc(3 * n * sizeof(double)); + // 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; @@ -177,9 +178,8 @@ static void fpli_free_cdllist(dlnode_t * head) { 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->r_next); - free(head[1].area); // free ->area and ->vol. + 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). } From 98129ef32900e968e84ebe256844c912ed8943d4 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 19:51:38 +0000 Subject: [PATCH 15/19] * hv.c (upper_bound): Move to hv_priv.h * hv4d_priv.h (update_bound_3d): Delete. Replaced by upper_bound(). --- c/hv.c | 8 -------- c/hv4d_priv.h | 13 +++---------- c/hv_priv.h | 10 ++++++++++ 3 files changed, 13 insertions(+), 18 deletions(-) diff --git a/c/hv.c b/c/hv.c index 2f36f0f80..df25be792 100644 --- a/c/hv.c +++ b/c/hv.c @@ -288,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, diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index cf4784694..16a810220 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -377,13 +377,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) -{ - for (int i = 0; i < 3; i++) - dest[i] = MAX(a[i], b[i]); -} - // FIXME: Move this function to hv.c #ifdef HV_RECURSIVE /** @@ -436,7 +429,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, return 0; } // x_aux is the coordinate-wise maximum between newpx and the_point_x. - update_bound_3d(x_aux, newpx, the_point_x); + upper_bound(x_aux, newpx, the_point_x, 3); newp_aux->x = x_aux; x_aux += 3; @@ -455,7 +448,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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. - update_bound_3d(x_aux, newpx, the_point_x); + upper_bound(x_aux, newpx, the_point_x, 3); x_aux += 3; c++; } @@ -518,7 +511,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, break; if (newp->ignore < 3) { // x_aux is the coordinate-wise maximum between newpx and the_point_x. - update_bound_3d(x_aux, newpx, 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)) { 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 From a91979c8952acbabe3816ff35bb6360e8efdea6c Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 22:22:37 +0000 Subject: [PATCH 16/19] Documentation. --- r/NEWS.md | 1 + 1 file changed, 1 insertion(+) 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 From 81b6ce1c33d6677ce8ddabe17e3ce3275fa6140f Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 8 Jan 2026 22:29:20 +0000 Subject: [PATCH 17/19] * c/hvc4d_priv.h: New file. --- c/Makefile | 3 +- c/hv.c | 2 +- c/hv4d_priv.h | 286 --------------------------------------------- c/hvc4d_priv.h | 311 +++++++++++++++++++++++++++++++++++++++++++++++++ c/libhv.mk | 2 +- 5 files changed, 315 insertions(+), 289 deletions(-) create mode 100644 c/hvc4d_priv.h 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/hv.c b/c/hv.c index df25be792..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 diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 16a810220..7c51c6c0b 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -171,91 +171,6 @@ 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) @@ -346,205 +261,4 @@ 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); -} - -// 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]); - - 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; - 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(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]; - } - } - - 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. -#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) { - // 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 // HV_RECURSIVE - #endif // _HV4D_PRIV_H diff --git a/c/hvc4d_priv.h b/c/hvc4d_priv.h new file mode 100644 index 000000000..4c3aa2d7e --- /dev/null +++ b/c/hvc4d_priv.h @@ -0,0 +1,311 @@ +#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" + +/** + 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; +} + +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); +} + +/** + 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; + 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(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]; + } + } + + 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. +#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) { + // 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 From 19b6389eb14ddd28f7ee2fa508880acce1bf1fa6 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 9 Jan 2026 09:49:10 +0000 Subject: [PATCH 18/19] * c/hvc4d_priv.h: Remove duplicated compare_lex_2d. --- c/hvc4d_priv.h | 13 +------------ c/sort.h | 5 ++--- 2 files changed, 3 insertions(+), 15 deletions(-) diff --git a/c/hvc4d_priv.h b/c/hvc4d_priv.h index 4c3aa2d7e..c30554d84 100644 --- a/c/hvc4d_priv.h +++ b/c/hvc4d_priv.h @@ -108,17 +108,6 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new return true; } -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) { @@ -130,7 +119,7 @@ lex_sort_equal_z_and_setup_nodes(dlnode_t * newp_aux, double * x_aux, size_t n) x += 3; } - qsort(scratch, n, sizeof(*scratch), compare_lex_2d); + qsort(scratch, n, sizeof(*scratch), cmp_ppdouble_asc_rev_2d); for (i = 0; i < n; i++) { newp_aux->x = scratch[i]; diff --git a/c/sort.h b/c/sort.h index ca211d5b4..224e0ee69 100644 --- a/c/sort.h +++ b/c/sort.h @@ -116,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 From 4bf704466354da7e51d5a4d38069fc0f130b9233 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 9 Jan 2026 13:35:42 +0000 Subject: [PATCH 19/19] * c/hvc4d_priv.h: Use lex_sort_buffer. --- c/hvc4d_priv.h | 73 +++++++++++++++++++++++++++++++------------------- c/sort.h | 2 +- 2 files changed, 46 insertions(+), 29 deletions(-) diff --git a/c/hvc4d_priv.h b/c/hvc4d_priv.h index c30554d84..08f5b793e 100644 --- a/c/hvc4d_priv.h +++ b/c/hvc4d_priv.h @@ -25,6 +25,44 @@ #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 @@ -108,28 +146,6 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new return true; } -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), cmp_ppdouble_asc_rev_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); -} - /** Compute the hv contribution of "the_point" in d=4 by iteratively computing the one contribution problem in d=3. */ @@ -195,6 +211,8 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, // 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) { @@ -211,7 +229,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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); + 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++; @@ -231,6 +249,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, } newp = newp->next[0]; } + free(sort_buffer.scratch); } dlnode_t * newp = the_point->next[1]; @@ -240,11 +259,9 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, dlnode_t * tp_prev_z = the_point->prev[0]; dlnode_t * tp_next_z = the_point->next[0]; // This call should 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 + _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; diff --git a/c/sort.h b/c/sort.h index 224e0ee69..9a59b2843 100644 --- a/c/sort.h +++ b/c/sort.h @@ -57,7 +57,7 @@ lexicographic_less_3d(const double * restrict a, const double * restrict b) 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])); + return a[1] < b[1] || (a[1] == b[1] && a[0] <= b[0]); } static inline bool