From 53ece9f87124f0c3cafc24b4c4c46e3020881824 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 2 Nov 2025 11:10:37 +0000 Subject: [PATCH 01/49] Start development of v0.2.0 --- Makefile | 2 +- c/Makefile | 2 +- python/pyproject.toml | 2 +- r/DESCRIPTION | 2 +- release.py | 2 ++ 5 files changed, 6 insertions(+), 4 deletions(-) diff --git a/Makefile b/Makefile index 89c2effa4..efb16736a 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -PACKAGEVERSION=0.1.9 +PACKAGEVERSION=0.2.0 .PHONY: default clean check test pre-commit diff --git a/c/Makefile b/c/Makefile index 68e1c08c5..b11a775ce 100644 --- a/c/Makefile +++ b/c/Makefile @@ -1,7 +1,7 @@ #----------------------------------------------------------------------- # Makefile for moocore #----------------------------------------------------------------------- -VERSION = 0.17.0 +VERSION = 0.17.0$(REVISION) S= ## Quiet / verbose output: diff --git a/python/pyproject.toml b/python/pyproject.toml index 6997bb6b1..56274773f 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -8,7 +8,7 @@ requires = [ [project] name = "moocore" -version = "0.1.9" +version = "0.2.0.dev0" description = "Core Algorithms for Multi-Objective Optimization" readme = "README.md" keywords = [ diff --git a/r/DESCRIPTION b/r/DESCRIPTION index 42d68c92c..d7f2ff64d 100644 --- a/r/DESCRIPTION +++ b/r/DESCRIPTION @@ -1,7 +1,7 @@ Package: moocore Type: Package Title: Core Mathematical Functions for Multi-Objective Optimization -Version: 0.1.9 +Version: 0.1.9.900 Authors@R: c(person("Manuel", "López-Ibáñez", role = c("aut", "cre"), email = "manuel.lopez-ibanez@manchester.ac.uk", comment = c(ORCID = "0000-0001-9974-1295")), diff --git a/release.py b/release.py index a2ea0cc24..3850724d0 100755 --- a/release.py +++ b/release.py @@ -111,6 +111,8 @@ def main(): update_package_version( "python/pyproject.toml", "version", replace=args.dev + ".dev0" ) + print("*** If happy with the changes, then do:") + print(f"git ci -a -m 'Start development of v{version}'") else: print(f"Current version is: {version}") print("Preparing for release") From 8f89f5e70ba151ffc10cd7f15e6672b0a26c3c07 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 31 Oct 2025 18:43:20 +0000 Subject: [PATCH 02/49] c/hv4d.c (restart_base_setup_z_and_closest): Simplify. Help autovectorization. Use _attr_optimize_finite_math. --- c/hv4d.c | 85 +++++++++++++++++++++++++++++++------------------------- 1 file changed, 47 insertions(+), 38 deletions(-) diff --git a/c/hv4d.c b/c/hv4d.c index 9029eec3e..87e0fe3ae 100644 --- a/c/hv4d.c +++ b/c/hv4d.c @@ -75,11 +75,11 @@ static void update_links(dlnode_t * restrict list, dlnode_t * restrict new) { assert(list+2 == list->prev[0]); - const double * newx = new->x; + const double newx[] = { new->x[0], new->x[1], new->x[2] }; dlnode_t * p = new->next[0]; - const double * px = p->x; - dlnode_t * stop = list+2; + const dlnode_t * stop = list+2; while (p != stop) { + const double * px = p->x; // px dominates newx (but not equal) if (px[0] <= newx[0] && px[1] <= newx[1] && (px[0] < newx[0] || px[1] < newx[1])) return; @@ -93,69 +93,78 @@ update_links(dlnode_t * restrict list, dlnode_t * restrict new) } else if (newx[0] < px[0] && lex_cmp_3d_102(newx, p->closest[1]->x)) { // newx[1] > px[1] p->closest[1] = new; } - } else if (newx[1] < px[1] && lex_cmp_3d_012(newx, p->closest[0]->x)) {//newx[0] > px[0] + } else if (newx[1] < px[1] && lex_cmp_3d_012(newx, p->closest[0]->x)) { //newx[0] > px[0] p->closest[0] = new; } p = p->next[0]; - px = p->x; } } // This does what setupZandClosest does while reconstructing L at z = new->x[2]. -__attribute__((hot)) static bool +_attr_optimize_finite_math +__attribute__((hot)) +static bool restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict new) { - const double * restrict newx = new->x; // FIXME: This is the most expensive function in the HV4D+ algorithm. + const double newx[] = { new->x[0], new->x[1], new->x[2], new->x[3] }; assert(list+1 == list->next[0]); - dlnode_t * closest1 = list; dlnode_t * closest0 = list+1; + dlnode_t * closest1 = list; const double * closest0x = closest0->x; const double * closest1x = closest1->x; - dlnode_t * p = list+1; - assert(p->next[0] == list->next[0]->next[0]); + dlnode_t * p = (list+1)->next[0]; + assert(p == list->next[0]->next[0]); restart_list_y(list); while (true) { - p = p->next[0]; const double * restrict px = p->x; // Help auto-vectorization. - bool p_leq_new_0 = px[0] <= newx[0]; - bool p_leq_new_1 = px[1] <= newx[1]; - bool p_leq_new_2 = px[2] <= newx[2]; - - if (p_leq_new_0 & p_leq_new_1 & p_leq_new_2) { + bool p_lt_new_0 = px[0] < newx[0]; + bool p_lt_new_1 = px[1] < newx[1]; + bool p_lt_new_2 = px[2] < newx[2]; + bool p_eq_new_0 = px[0] == newx[0]; + bool p_eq_new_1 = px[1] == newx[1]; + bool p_eq_new_2 = px[2] == newx[2]; + bool p_leq_new_0 = p_lt_new_0 | p_eq_new_0; + bool p_leq_new_1 = p_lt_new_1 | p_eq_new_1; + bool p_leq_new_2 = p_lt_new_2 | p_eq_new_2; + + // if (weakly_dominates(px, newx, 3)) { // Slower + if (unlikely(p_leq_new_0 & p_leq_new_1 & p_leq_new_2)) { //new->ndomr++; assert(weakly_dominates(px, newx, 4)); return false; } - if (!lexicographic_less_3d(px, newx)) - break; + //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))))) { + new->closest[0] = closest0; + new->closest[1] = closest1; + new->prev[0] = p->prev[0]; + new->next[0] = p; + return true; + } - // reconstruct + // FIXME: Can we move reconstruct() after setup_z_and_closest() ? + // reconstruct() set_cnext_to_closest(p); p->cnext[0]->cnext[1] = p; p->cnext[1]->cnext[0] = p; - // setup_z_and_closest - assert(px[0] > newx[0] || px[1] > newx[1]); - if (px[1] < newx[1] && (px[0] < closest0x[0] || (px[0] == closest0x[0] && px[1] < closest0x[1]))) { + // setup_z_and_closest() + ASSUME(!p_leq_new_0 || !p_leq_new_1); + if (p_lt_new_1 && (px[0] < closest0x[0] || (px[0] == closest0x[0] && px[1] < closest0x[1]))) { closest0 = p; closest0x = px; - } else if (px[0] < newx[0] && (px[1] < closest1x[1] || (px[1] == closest1x[1] && px[0] < closest1x[0]))) { + } else if (p_lt_new_0 && (px[1] < closest1x[1] || (px[1] == closest1x[1] && px[0] < closest1x[0]))) { closest1 = p; closest1x = px; } + p = p->next[0]; } - - new->closest[0] = closest0; - new->closest[1] = closest1; - - new->prev[0] = p->prev[0]; - new->next[0] = p; - return true; + unreachable(); } // FIXME: This is very similar to the loop in hvc3d_list() but it doesn't use p->last_slice_z @@ -167,19 +176,17 @@ one_contribution_3d(dlnode_t * restrict new) // if newx[0] == new->cnext[0]->x[0], the first area is zero double area = compute_area_no_inners(newx, new->cnext[0], 1); double volume = 0; - dlnode_t * p = new; - const double * px = p->x; - assert(!weakly_dominates(p->next[0]->x, newx, 4)); + double lastz = newx[2]; + dlnode_t * p = new->next[0]; + assert(!weakly_dominates(p->x, newx, 4)); while (true) { - double lastz = px[2]; - p = p->next[0]; - px = p->x; + const double * px = p->x; volume += area * (px[2] - lastz); if (px[0] <= newx[0] && px[1] <= newx[1]) return volume; - assert(px[0] > newx[0] || px[1] > newx[1]); + ASSUME(px[0] > newx[0] || px[1] > newx[1]); assert(!weakly_dominates(px, p->next[0]->x, 4)); set_cnext_to_closest(p); @@ -209,8 +216,10 @@ one_contribution_3d(dlnode_t * restrict new) p->cnext[1]->cnext[0] = p; p->cnext[0]->cnext[1] = p; } + lastz = px[2]; + p = p->next[0]; } - return volume; + unreachable(); } /* Compute the hypervolume indicator in d=4 by iteratively computing the one From 09b418656e7d7c12f4e61fa4935a44f3301509d9 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 1 Nov 2025 22:20:19 +0000 Subject: [PATCH 03/49] c/eaf.c (eaf_check_polygons): Do not call polygons_intersect() twice. --- c/eaf.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/c/eaf.c b/c/eaf.c index 16e492b5d..a924ae4b5 100644 --- a/c/eaf.c +++ b/c/eaf.c @@ -563,7 +563,7 @@ eaf_check_polygons(eaf_polygon_t *p, int nobj) polygon_print(pi, nobj); polygon_print(pj, nobj); #endif - assert(!polygons_intersect(pi, pj, nobj)); + assert(false); } pj = next_polygon(pj, nobj, end); } From faa5579da3e0c1f165a4cefd5c0a81c5ab9d6a2e Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 2 Nov 2025 10:54:03 +0000 Subject: [PATCH 04/49] README.md, c/README.md: Fix testsuite badge. --- README.md | 2 +- c/README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 0eda9d0be..25644ab3e 100644 --- a/README.md +++ b/README.md @@ -79,5 +79,5 @@ The following projects currently use `moocore`: [r-universe-build-link]: https://github.com/r-universe/multi-objective/actions/workflows/build.yml [r-universe-version]: https://multi-objective.r-universe.dev/moocore [r-universe-version-badge]: https://multi-objective.r-universe.dev/badges/moocore -[testsuite-badge]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml/badge.svg?event=push +[testsuite-badge]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml/badge.svg [testsuite-link]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml diff --git a/c/README.md b/c/README.md index 85d49e435..6e5879ff7 100644 --- a/c/README.md +++ b/c/README.md @@ -348,5 +348,5 @@ See the [LICENSE](/LICENSE) and [COPYRIGHTS](/r/inst/COPYRIGHTS) files. [r-moocore-cran]: https://cran.r-project.org/package=moocore [r-moocore-github]: https://github.com/multi-objective/moocore/tree/main/r#readme [r-moocore-homepage]: https://multi-objective.github.io/moocore/r/ -[testsuite-badge]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml/badge.svg?event=push +[testsuite-badge]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml/badge.svg [testsuite-link]: https://github.com/multi-objective/testsuite/actions/workflows/moocore.yml From 2b29d737fccead4588a1c55f9d8baf6c6aa82900 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 2 Nov 2025 12:28:28 +0000 Subject: [PATCH 05/49] * python/src/moocore/_utils.py (_get_seed_for_c): New. * python/src/moocore/_moocore.py: Use it. --- python/src/moocore/_moocore.py | 12 +++--------- python/src/moocore/_utils.py | 6 ++++++ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index 578be70bc..c25d146de 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -23,6 +23,7 @@ np1d_to_int_array, atleast_1d_of_length_n, is_integer_value, + _get_seed_for_c, ) ## The CFFI library is used to create C bindings. @@ -938,11 +939,7 @@ def hv_approx( raise ValueError(f"nsamples must be an integer value: {nsamples}") nsamples = ffi.cast("uint_fast32_t", nsamples) if method == "DZ2019-MC": - if not is_integer_value(seed): - seed = np.random.default_rng(seed).integers( - 2**32 - 2, dtype=np.uint32 - ) - seed = ffi.cast("uint32_t", seed) + seed = _get_seed_for_c(seed) hv = lib.hv_approx_normal( data_p, nobj, npoints, ref, maximise, nsamples, seed ) @@ -2545,14 +2542,11 @@ def whv_hype( ref[maximise] = -ref[maximise] ideal[maximise] = -ideal[maximise] - if not is_integer_value(seed): - seed = np.random.default_rng(seed).integers(2**32 - 2, dtype=np.uint32) - data_p, npoints, nobj = np2d_to_double_array(data) ref = ffi.from_buffer("double []", ref) ideal = ffi.from_buffer("double []", ideal) # FIXME: Check ranges. - seed = ffi.cast("uint32_t", seed) + seed = _get_seed_for_c(seed) nsamples = ffi.cast("int", nsamples) if dist == "uniform": diff --git a/python/src/moocore/_utils.py b/python/src/moocore/_utils.py index 26268e744..f7f403984 100644 --- a/python/src/moocore/_utils.py +++ b/python/src/moocore/_utils.py @@ -71,3 +71,9 @@ def is_integer_value(n): return False except TypeError: return False + + +def _get_seed_for_c(seed): + if not is_integer_value(seed): + seed = np.random.default_rng(seed).integers(2**32 - 2, dtype=np.uint32) + return ffi.cast("uint32_t", seed) From 2224a59384641245633851f4583fad0a367d3ebe Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 2 Nov 2025 12:29:13 +0000 Subject: [PATCH 06/49] * python/tests/test_moocore.py (test_hv_approx_default_seed): New test. (check_hvc): check also maximise=True. --- python/tests/test_moocore.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/python/tests/test_moocore.py b/python/tests/test_moocore.py index 6e8463536..92d5f634f 100644 --- a/python/tests/test_moocore.py +++ b/python/tests/test_moocore.py @@ -372,19 +372,33 @@ def test_hv_approx(dim): np.testing.assert_approx_equal(true_hv, appr_hv, significant=signif) +def test_hv_approx_default_seed(): + x = np.full((1, 5), 0.5) + ref = np.full(5, 0.0) + true_hv = moocore.hypervolume(x, ref=ref, maximise=True) + appr_hv = moocore.hv_approx( + x, ref=ref, method="DZ2019-MC", seed=None, maximise=True + ) + ## Precision goes down significantly with higher dimensions. + signif = 2 + np.testing.assert_approx_equal(true_hv, appr_hv, significant=signif) + + def check_hvc(points, ref, err_msg): hvc = moocore.hv_contributions(points, ref=ref) - nondom = moocore.is_nondominated(points, keep_weakly=True) - points = points[nondom, :] - hv_total = moocore.hypervolume(points, ref=ref) - true_hvc = np.zeros(len(nondom)) - true_hvc[nondom] = hv_total - np.array( + is_nondom = moocore.is_nondominated(points, keep_weakly=True) + nondom = points[is_nondom, :] + hv_total = moocore.hypervolume(nondom, ref=ref) + true_hvc = np.zeros(len(is_nondom)) + true_hvc[is_nondom] = hv_total - np.array( [ - moocore.hypervolume(np.delete(points, i, axis=0), ref=ref) - for i in range(len(points)) + moocore.hypervolume(np.delete(nondom, i, axis=0), ref=ref) + for i in range(len(nondom)) ] ) assert_allclose(true_hvc, hvc, err_msg=err_msg) + hvc = moocore.hv_contributions(-points, ref=-ref, maximise=True) + assert_allclose(true_hvc, hvc, err_msg=err_msg) @pytest.mark.parametrize("dim", range(2, 5)) From dd1ddbb5581a15bac44e16432eefc606dff1e5bf Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sun, 2 Nov 2025 14:56:29 +0000 Subject: [PATCH 07/49] * c/hv_priv.h (init_sentinels): Make sure HV_DIMENSION is correct. --- c/hv_priv.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/c/hv_priv.h b/c/hv_priv.h index 6dadac436..aa6569314 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -145,11 +145,13 @@ init_sentinels(dlnode_t * restrict list, const double * restrict ref) #if HV_DIMENSION == 3 -DBL_MAX, ref[1], -DBL_MAX, // Sentinel 1 ref[0], -DBL_MAX, -DBL_MAX, // Sentinel 2 - -DBL_MAX, -DBL_MAX, ref[2] // Sentinel 2 -#else + -DBL_MAX, -DBL_MAX, ref[2] // Sentinel 3 +#elif HV_DIMENSION == 4 -DBL_MAX, ref[1], -DBL_MAX, -DBL_MAX, // Sentinel 1 ref[0], -DBL_MAX, -DBL_MAX, -DBL_MAX, // Sentinel 2 - -DBL_MAX, -DBL_MAX, ref[2], ref[3] // Sentinel 2 + -DBL_MAX, -DBL_MAX, ref[2], ref[3] // Sentinel 3 +#else +# error "Unknown HV_DIMENSION" #endif }; From 9f4c5de7c0174e84e09018d4fba691283ccc4ae7 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 3 Nov 2025 09:57:33 +0000 Subject: [PATCH 08/49] * c/mt19937/mt19937.c (mt19937_init_by_array,init_genrand): Remove unused. * c/mt19937/mt19937.h: Remove declarations. --- c/mt19937/mt19937.c | 98 --------------------------------------------- c/mt19937/mt19937.h | 3 -- 2 files changed, 101 deletions(-) diff --git a/c/mt19937/mt19937.c b/c/mt19937/mt19937.c index 0ca524ffc..8da3db733 100644 --- a/c/mt19937/mt19937.c +++ b/c/mt19937/mt19937.c @@ -15,104 +15,6 @@ void mt19937_seed(mt19937_state *state, uint32_t seed) { state->pos = RK_STATE_LEN; } -/* initializes mt[RK_STATE_LEN] with a seed */ -static void init_genrand(mt19937_state *state, uint32_t s) { - int mti; - uint32_t *mt = state->key; - - mt[0] = s & 0xffffffffUL; - for (mti = 1; mti < RK_STATE_LEN; mti++) { - /* - * See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. - * In the previous versions, MSBs of the seed affect - * only MSBs of the array mt[]. - * 2002/01/09 modified by Makoto Matsumoto - */ - /* This code triggers conversions between uint32_t and unsigned long, - which various compilers warn about. */ -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wconversion" -#if defined(__clang__) -# pragma clang diagnostic push -# pragma clang diagnostic ignored "-Wimplicit-int-conversion" -#endif - mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti); -# pragma GCC diagnostic pop -#if defined(__clang__) -# pragma clang diagnostic pop -#endif - /* for > 32 bit machines */ - mt[mti] &= 0xffffffffUL; - } - state->pos = mti; - return; -} - -/* - * initialize by an array with array-length - * init_key is the array for initializing keys - * key_length is its length - */ -void mt19937_init_by_array(mt19937_state *state, uint32_t *init_key, - int key_length) { - /* was signed in the original code. RDH 12/16/2002 */ - int i = 1; - int j = 0; - uint32_t *mt = state->key; - int k; - - init_genrand(state, 19650218UL); - k = (RK_STATE_LEN > key_length ? RK_STATE_LEN : key_length); - for (; k; k--) { -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wconversion" -#if defined(__clang__) -# pragma clang diagnostic push -# pragma clang diagnostic ignored "-Wimplicit-int-conversion" -#endif - /* non linear */ - mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1664525UL)) + init_key[j] + j; -# pragma GCC diagnostic pop -#if defined(__clang__) -# pragma clang diagnostic pop -#endif - /* for > 32 bit machines */ - mt[i] &= 0xffffffffUL; - i++; - j++; - if (i >= RK_STATE_LEN) { - mt[0] = mt[RK_STATE_LEN - 1]; - i = 1; - } - if (j >= key_length) { - j = 0; - } - } - for (k = RK_STATE_LEN - 1; k; k--) { - /* This code triggers conversions between uint32_t and unsigned long, - which various compilers warn about. */ -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wconversion" -#if defined(__clang__) -# pragma clang diagnostic push -# pragma clang diagnostic ignored "-Wimplicit-int-conversion" -#endif - mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * 1566083941UL)) - i; /* non linear */ -# pragma GCC diagnostic pop -#if defined(__clang__) -# pragma clang diagnostic pop -#endif - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - if (i >= RK_STATE_LEN) { - mt[0] = mt[RK_STATE_LEN - 1]; - i = 1; - } - } - - mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ -} - void mt19937_gen(mt19937_state *state) { uint32_t y; int i; diff --git a/c/mt19937/mt19937.h b/c/mt19937/mt19937.h index 368fd799a..16c09c230 100644 --- a/c/mt19937/mt19937.h +++ b/c/mt19937/mt19937.h @@ -43,9 +43,6 @@ static inline uint32_t mt19937_next(mt19937_state *state) { return y; } -extern void mt19937_init_by_array(mt19937_state *state, uint32_t *init_key, - int key_length); - static inline uint64_t mt19937_next64(mt19937_state *state) { return (uint64_t)mt19937_next(state) << 32 | mt19937_next(state); } From 5d5a7df8c473f6b7fa22a693ceb2251fe4057079 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 3 Nov 2025 17:47:16 +0000 Subject: [PATCH 09/49] c/Makefile (clean): Clean .i and .s. --- c/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/c/Makefile b/c/Makefile index b11a775ce..e8dab4a06 100644 --- a/c/Makefile +++ b/c/Makefile @@ -224,6 +224,7 @@ $(CXXSHLIB): $(LIBHV_OBJS) $(CXX) $(CFLAGS) $(CPPFLAGS) $(LDFLAGS) -o $@ $^ $(SHLIB_LDFLAGS) clean: + @-$(RM) *.s *.i @-$(RM) config.status config.log @-$(RM) Hypervolume_MEX.mex @-$(RM) $(SHLIB) $(CXXSHLIB) From 02326474cdc33573c40621c26b7d15a281db1b6e Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 3 Nov 2025 13:20:08 +0000 Subject: [PATCH 10/49] c/epsilon.c: Fix error message. --- c/epsilon.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/c/epsilon.c b/c/epsilon.c index eb2b742a6..e0880a2e5 100644 --- a/c/epsilon.c +++ b/c/epsilon.c @@ -143,7 +143,7 @@ do_file (const char *filename, double *reference, size_t reference_size, // time_elapsed = Timer_elapsed_virtual (); fprintf (outfile, indicator_printf_format "\n", epsilon); if ((additive_flag && epsilon < 0) || (!additive_flag && epsilon < 1)) { - errprintf ("%s: some points are not dominated by the reference set", + errprintf ("%s: some points are not dominated by the reference set", filename); exit (EXIT_FAILURE); }/* From b899ffcc8f149ecfff17f67220d1e81cb270aa39 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 3 Nov 2025 13:26:51 +0000 Subject: [PATCH 11/49] c/epsilon.h: Use _attr_optimize_finite_math. Improve function specialization and vectorization. --- c/epsilon.h | 131 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 93 insertions(+), 38 deletions(-) diff --git a/c/epsilon.h b/c/epsilon.h index 552da5be1..49a716677 100644 --- a/c/epsilon.h +++ b/c/epsilon.h @@ -14,11 +14,12 @@ static inline bool all_positive(const double * restrict points, size_t size, dimension_t dim) { - ASSUME(dim <= 32); - for (size_t a = 0; a < size; a++) - for (dimension_t d = 0; d < dim; d++) - if (points[a * dim + d] <= 0) - return false; + ASSUME(size > 0); + ASSUME(1 <= dim && dim <= 32); + const size_t len = size * dim; + for (size_t a = 0; a < len; a++) + if (points[a] <= 0) + return false; return true; } @@ -43,19 +44,19 @@ all_positive(const double * restrict points, size_t size, dimension_t dim) for negative values doesn't make sense. */ -static inline double -eps_value_(bool do_ratio, double a, double b) -{ - return do_ratio ? a / b : a - b; -} +_attr_optimize_finite_math static inline double epsilon_helper_(bool do_mult, const enum objs_agree_t agree, const signed char * restrict minmax, dimension_t dim, const double * restrict points_a, size_t size_a, const double * restrict points_b, size_t size_b) { +// Converting this macro to an inline function hinders vectorization. +#define eps_value_(X,Y) (do_mult ? ((X) / (Y)) : ((X) - (Y))) + ASSUME(2 <= dim && dim <= 32); + ASSUME(size_a > 0 && size_b > 0); ASSUME(agree == AGREE_MINIMISE || agree == AGREE_MAXIMISE || agree == AGREE_NONE); ASSUME(minmax == NULL || agree != AGREE_NONE); double epsilon = do_mult ? 0 : -INFINITY; @@ -65,29 +66,25 @@ epsilon_helper_(bool do_mult, const enum objs_agree_t agree, const double * restrict pb = &points_b[b * dim]; for (size_t a = 0; a < size_a; a++) { const double * restrict pa = &points_a[a * dim]; - double epsilon_max; - if (agree == AGREE_NONE) { - epsilon_max = MAX(minmax[0] * eps_value_(do_mult, pb[0], pa[0]), - minmax[1] * eps_value_(do_mult, pb[1], pa[1])); - if (epsilon_max >= epsilon_min) - continue; - for (dimension_t d = 2; d < dim; d++) { - double epsilon_temp = minmax[d] * eps_value_(do_mult, pb[d], pa[d]); - epsilon_max = MAX(epsilon_max, epsilon_temp); - } - } else { - epsilon_max = (agree == AGREE_MINIMISE) - ? MAX(eps_value_(do_mult, pa[0], pb[0]), eps_value_(do_mult, pa[1], pb[1])) - : MAX(eps_value_(do_mult, pb[0], pa[0]), eps_value_(do_mult, pb[1], pa[1])); - if (epsilon_max >= epsilon_min) - continue; - for (dimension_t d = 2; d < dim; d++) { - double epsilon_temp = (agree == AGREE_MINIMISE) - ? eps_value_(do_mult, pa[d], pb[d]) - : eps_value_(do_mult, pb[d], pa[d]); - epsilon_max = MAX(epsilon_max, epsilon_temp); - } + double epsilon_max = (agree == AGREE_NONE) + ? MAX(minmax[0] * eps_value_(pb[0], pa[0]), + minmax[1] * eps_value_(pb[1], pa[1])) + : ((agree == AGREE_MINIMISE) + ? MAX(eps_value_(pa[0], pb[0]), eps_value_(pa[1], pb[1])) + : MAX(eps_value_(pb[0], pa[0]), eps_value_(pb[1], pa[1]))); + + if (epsilon_max >= epsilon_min) + continue; + + for (dimension_t d = 2; d < dim; d++) { + double epsilon_temp = (agree == AGREE_NONE) + ? minmax[d] * eps_value_(pb[d], pa[d]) + : ((agree == AGREE_MINIMISE) + ? eps_value_(pa[d], pb[d]) + : eps_value_(pb[d], pa[d])); + epsilon_max = MAX(epsilon_max, epsilon_temp); } + if (epsilon_max <= epsilon) { skip_max = true; break; @@ -100,6 +97,63 @@ epsilon_helper_(bool do_mult, const enum objs_agree_t agree, return epsilon; } +_attr_optimize_finite_math +static inline double +epsilon_mult_agree_none(const signed char * restrict minmax, dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/true, AGREE_NONE, minmax, dim, points_a, size_a, points_b, size_b); +} + +_attr_optimize_finite_math +static inline double +epsilon_mult_agree_min(dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/true, AGREE_MINIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); +} + +_attr_optimize_finite_math +static inline double +epsilon_mult_agree_max(dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/true, AGREE_MAXIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); +} + + +_attr_optimize_finite_math +static inline double +epsilon_addi_agree_none(const signed char * restrict minmax, dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/false, AGREE_NONE, minmax, dim, points_a, size_a, points_b, size_b); +} + +_attr_optimize_finite_math +static inline double +epsilon_addi_agree_min(dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/false, AGREE_MINIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); +} + +_attr_optimize_finite_math +static inline double +epsilon_addi_agree_max(dimension_t dim, + const double * restrict points_a, size_t size_a, + const double * restrict points_b, size_t size_b) +{ + return epsilon_helper_(/* do_mult=*/false, AGREE_MAXIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); +} + + +_attr_optimize_finite_math static inline double epsilon_mult_minmax(const signed char * restrict minmax, dimension_t dim, const double * restrict points_a, size_t size_a, @@ -114,14 +168,15 @@ epsilon_mult_minmax(const signed char * restrict minmax, dimension_t dim, // This forces the compiler to generate three specialized versions of the function. switch (check_all_minimize_maximize(minmax, dim)) { case AGREE_MINIMISE: - return epsilon_helper_(/* do_mult=*/true, AGREE_MINIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); + return epsilon_mult_agree_min(dim, points_a, size_a, points_b, size_b); case AGREE_MAXIMISE: - return epsilon_helper_(/* do_mult=*/true, AGREE_MAXIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); + return epsilon_mult_agree_max(dim, points_a, size_a, points_b, size_b); default: - return epsilon_helper_(/* do_mult=*/true, AGREE_NONE, minmax, dim, points_a, size_a, points_b, size_b); + return epsilon_mult_agree_none(minmax, dim, points_a, size_a, points_b, size_b); } } +_attr_optimize_finite_math static inline double epsilon_additive_minmax(const signed char * restrict minmax, dimension_t dim, const double * restrict points_a, size_t size_a, @@ -130,11 +185,11 @@ epsilon_additive_minmax(const signed char * restrict minmax, dimension_t dim, // This forces the compiler to generate three specialized versions of the function. switch (check_all_minimize_maximize(minmax, dim)) { case AGREE_MINIMISE: - return epsilon_helper_(/* do_mult=*/false, AGREE_MINIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); + return epsilon_addi_agree_min(dim, points_a, size_a, points_b, size_b); case AGREE_MAXIMISE: - return epsilon_helper_(/* do_mult=*/false, AGREE_MAXIMISE, /*minmax=*/NULL, dim, points_a, size_a, points_b, size_b); + return epsilon_addi_agree_max(dim, points_a, size_a, points_b, size_b); default: - return epsilon_helper_(/* do_mult=*/false, AGREE_NONE, minmax, dim, points_a, size_a, points_b, size_b); + return epsilon_addi_agree_none(minmax, dim, points_a, size_a, points_b, size_b); } } From 5ee51d1bcbbec049ebbce852a8c3949e7520c514 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 3 Nov 2025 16:58:52 +0000 Subject: [PATCH 12/49] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/lorenzwalthert/precommit: v0.4.3.9014 → v0.4.3.9017](https://github.com/lorenzwalthert/precommit/compare/v0.4.3.9014...v0.4.3.9017) - [github.com/tox-dev/tox-ini-fmt: 1.6.0 → 1.7.0](https://github.com/tox-dev/tox-ini-fmt/compare/1.6.0...1.7.0) - [github.com/tox-dev/pyproject-fmt: v2.6.0 → v2.11.0](https://github.com/tox-dev/pyproject-fmt/compare/v2.6.0...v2.11.0) - [github.com/astral-sh/ruff-pre-commit: v0.13.0 → v0.14.3](https://github.com/astral-sh/ruff-pre-commit/compare/v0.13.0...v0.14.3) - [github.com/sphinx-contrib/sphinx-lint: v1.0.0 → v1.0.1](https://github.com/sphinx-contrib/sphinx-lint/compare/v1.0.0...v1.0.1) --- .pre-commit-config.yaml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 9f2d2bd4d..7e6350ad5 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ ci: repos: - repo: https://github.com/lorenzwalthert/precommit - rev: v0.4.3.9014 + rev: v0.4.3.9017 hooks: - id: parsable-R - id: no-browser-statement @@ -43,19 +43,19 @@ repos: exclude: '(\.Rd|python/doc/source/reference/.*|test-doctest-.*)' - repo: https://github.com/tox-dev/tox-ini-fmt - rev: "1.6.0" + rev: "1.7.0" hooks: - id: tox-ini-fmt - repo: https://github.com/tox-dev/pyproject-fmt - rev: "v2.6.0" + rev: "v2.11.0" hooks: - id: pyproject-fmt additional_dependencies: ["tox>=4.12.1"] - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.13.0 + rev: v0.14.3 hooks: # Run the formatter. - id: ruff-format @@ -67,6 +67,6 @@ repos: require_serial: true - repo: https://github.com/sphinx-contrib/sphinx-lint - rev: v1.0.0 + rev: v1.0.1 hooks: - id: sphinx-lint From ccb2dd825402444a0ee40d4ea5395e3f14a6d580 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 3 Nov 2025 17:00:13 +0000 Subject: [PATCH 13/49] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- python/pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/python/pyproject.toml b/python/pyproject.toml index 56274773f..da1d4b9f0 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -33,6 +33,7 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", ] dependencies = [ "cffi>=1.17.1", From d443ecf0b6cccaaad36dfe420cc967b4c01695e1 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 4 Nov 2025 09:33:14 +0000 Subject: [PATCH 14/49] r/DESCRIPTION: Edit Description. --- c/libmoocore-config.h | 2 +- r/DESCRIPTION | 2 +- r/man/moocore-package.Rd | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/c/libmoocore-config.h b/c/libmoocore-config.h index 22a9a495c..c64b3c86c 100644 --- a/c/libmoocore-config.h +++ b/c/libmoocore-config.h @@ -2,7 +2,7 @@ # define LIBMOOCORE_CONFIG_H_ #ifndef R_PACKAGE -# if defined(_WIN32) && !(defined(__MINGW__) || defined(__MINGW32__) || defined(__MINGW64__)) +# if (defined(_WIN32) || defined(__CYGWIN__)) && !(defined(__MINGW__) || defined(__MINGW32__) || defined(__MINGW64__)) # ifndef MOOCORE_STATIC_LIB # ifdef MOOCORE_SHARED_LIB # define MOOCORE_API extern __declspec(dllexport) diff --git a/r/DESCRIPTION b/r/DESCRIPTION index d7f2ff64d..982a5835b 100644 --- a/r/DESCRIPTION +++ b/r/DESCRIPTION @@ -15,7 +15,7 @@ Authors@R: c(person("Manuel", "López-Ibáñez", role = c("aut", "cre"), person("Jean-Sebastien", "Roy", role = "cph", comment = "mt19937 library"), person("Makoto", "Matsumoto", role = "cph", comment = "mt19937 library"), person("Takuji", "Nishimura", role = "cph", comment = "mt19937 library")) -Description: Fast implementation of mathematical operations and performance metrics for multi-objective optimization, including filtering and ranking of dominated vectors according to Pareto optimality, computation of the empirical attainment function, V.G. da Fonseca, C.M. Fonseca, A.O. Hall (2001) , hypervolume metric, C.M. Fonseca, L. Paquete, M. López-Ibáñez (2006) , epsilon indicator, inverted generational distance, and Vorob'ev threshold, expectation and deviation, M. Binois, D. Ginsbourger, O. Roustant (2015) , among others. +Description: Fast implementations of mathematical operations and performance metrics for multi-objective optimization, including filtering and ranking of dominated vectors according to Pareto optimality, hypervolume metric, C.M. Fonseca, L. Paquete, M. López-Ibáñez (2006) , epsilon indicator, inverted generational distance, computation of the empirical attainment function, V.G. da Fonseca, C.M. Fonseca, A.O. Hall (2001) , and Vorob'ev threshold, expectation and deviation, M. Binois, D. Ginsbourger, O. Roustant (2015) , among others. Depends: R (>= 4.0) Imports: matrixStats, diff --git a/r/man/moocore-package.Rd b/r/man/moocore-package.Rd index fa915ec5c..417dca2c0 100644 --- a/r/man/moocore-package.Rd +++ b/r/man/moocore-package.Rd @@ -6,7 +6,7 @@ \alias{moocore-package} \title{moocore: Core Mathematical Functions for Multi-Objective Optimization} \description{ -Fast implementation of mathematical operations and performance metrics for multi-objective optimization, including filtering and ranking of dominated vectors according to Pareto optimality, computation of the empirical attainment function, V.G. da Fonseca, C.M. Fonseca, A.O. Hall (2001) \doi{10.1007/3-540-44719-9_15}, hypervolume metric, C.M. Fonseca, L. Paquete, M. López-Ibáñez (2006) \doi{10.1109/CEC.2006.1688440}, epsilon indicator, inverted generational distance, and Vorob'ev threshold, expectation and deviation, M. Binois, D. Ginsbourger, O. Roustant (2015) \doi{10.1016/j.ejor.2014.07.032}, among others. +Fast implementations of mathematical operations and performance metrics for multi-objective optimization, including filtering and ranking of dominated vectors according to Pareto optimality, hypervolume metric, C.M. Fonseca, L. Paquete, M. López-Ibáñez (2006) \doi{10.1109/CEC.2006.1688440}, epsilon indicator, inverted generational distance, computation of the empirical attainment function, V.G. da Fonseca, C.M. Fonseca, A.O. Hall (2001) \doi{10.1007/3-540-44719-9_15}, and Vorob'ev threshold, expectation and deviation, M. Binois, D. Ginsbourger, O. Roustant (2015) \doi{10.1016/j.ejor.2014.07.032}, among others. } \seealso{ Useful links: From 8290a6656b2f6753799d9390355d3ac02faa6703 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 4 Nov 2025 12:13:27 +0000 Subject: [PATCH 15/49] c/hvapprox.c: Simplify ASSUME conditions. --- c/hvapprox.c | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/c/hvapprox.c b/c/hvapprox.c index 446aa687c..b69cd1013 100644 --- a/c/hvapprox.c +++ b/c/hvapprox.c @@ -29,7 +29,7 @@ transform_and_filter(const double * restrict data, size_t * restrict npoints_p, const bool * restrict maximise) { size_t npoints = *npoints_p; - double * points = malloc(dim * npoints * sizeof(double)); + double * points = malloc(dim * npoints * sizeof(*points)); size_t i, j; // Transform points (ref - points) for (i = 0, j = 0; i < npoints; i++) { @@ -54,14 +54,14 @@ transform_and_filter(const double * restrict data, size_t * restrict npoints_p, } _attr_optimize_finite_math // Required so that GCC will vectorize the inner loop. -static inline double +static double get_expected_value(const double * restrict points, size_t npoints, dimension_t dim, const double * restrict w) { ASSUME(1 <= dim && dim <= 32); ASSUME(npoints > 0); // points >= 0 && w >=0 so max_s_w cannot be < 0. - double max_s_w = 0; + double max_s_w = -INFINITY; for (size_t i = 0; i < npoints; i++) { const double * restrict p = points + i * dim; double min_ratio = p[0] * w[0]; @@ -71,6 +71,7 @@ get_expected_value(const double * restrict points, size_t npoints, } max_s_w = MAX(max_s_w, min_ratio); } + ASSUME(max_s_w >= 0); return pow_uint(max_s_w, dim); } @@ -160,7 +161,7 @@ static const long double sphere_area_div_2_pow_d_times_d[] = { _attr_pure_func static double euclidean_norm(const double * restrict w, dimension_t dim) { - ASSUME(dim >= 2 && dim <= 32); + ASSUME(2 <= dim && dim <= 32); double norm = (w[0] * w[0]) + (w[1] * w[1]); for (dimension_t k = 2; k < dim; k++) norm += w[k] * w[k]; @@ -530,8 +531,6 @@ solve_inverse_int_of_power_sin(long double theta, dimension_t dim) static long double * compute_int_all(dimension_t dm1) { - ASSUME(dm1 >= 1); - ASSUME(dm1 < 32); long double * int_all = malloc(dm1 * sizeof(long double)); dimension_t i; DEBUG2_PRINT("int_all[%u] =", (unsigned int)dm1); @@ -560,8 +559,7 @@ static void compute_theta(long double * restrict theta, dimension_t dim, const long double * restrict int_all) { - ASSUME(dim >= 2); - ASSUME(dim <= 32); + ASSUME(2 <= dim && dim <= 32); for (dimension_t j = 0; j < dim - 1; j++) { // We multiply here because we computed 1 / int_all[j] before. theta[j] = solve_inverse_int_of_power_sin(theta[j] * int_all[(dim - 2) - j], @@ -573,8 +571,7 @@ static void compute_hua_wang_direction(double * restrict direction, dimension_t dim, const long double * restrict theta) { - ASSUME(dim >= 2); - ASSUME(dim <= 32); + ASSUME(2 <= dim && dim <= 32); dimension_t k, j; direction[0] = STATIC_CAST(double, sinl(theta[0])); for (k = 1; k < dim - 1; k++) From f28adef8841d49e95e6106316f9df1468379350d Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 4 Nov 2025 14:50:12 +0000 Subject: [PATCH 16/49] README.md: Mention QuantumAnnealingMOO. --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 25644ab3e..cd5b33bae 100644 --- a/README.md +++ b/README.md @@ -53,6 +53,8 @@ The following projects currently use `moocore`: - [`pymoo`](https://pymoo.org/) - [`Liger`](https://www.ide.uk/project-liger) - [`DESDEO`](https://github.com/industrial-optimization-group/DESDEO/) + - [QuantumAnnealingMOO](https://github.com/andrew-d-king/QuantumAnnealingMOO) (D-Wave Quantum Inc) + [c-build-badge]: https://github.com/multi-objective/moocore/actions/workflows/C.yml/badge.svg?event=push From 958a0909d883a71a363a1112676f3670b866056c Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 5 Nov 2025 19:13:12 +0000 Subject: [PATCH 17/49] README.md: Mention IOHinspector --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index cd5b33bae..eb84f8743 100644 --- a/README.md +++ b/README.md @@ -53,6 +53,7 @@ The following projects currently use `moocore`: - [`pymoo`](https://pymoo.org/) - [`Liger`](https://www.ide.uk/project-liger) - [`DESDEO`](https://github.com/industrial-optimization-group/DESDEO/) + - [`IOHInspector`](https://github.com/IOHprofiler/IOHinspector?tab=readme-ov-file) - [QuantumAnnealingMOO](https://github.com/andrew-d-king/QuantumAnnealingMOO) (D-Wave Quantum Inc) From 3dff60b279199bf0376634cdd7806ce9db4a355d Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 6 Nov 2025 09:25:54 +0000 Subject: [PATCH 18/49] * python/doc/source/whatsnew/index.rst: Fix link. --- python/doc/source/whatsnew/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/doc/source/whatsnew/index.rst b/python/doc/source/whatsnew/index.rst index 293ea7de0..9dd8b7b32 100644 --- a/python/doc/source/whatsnew/index.rst +++ b/python/doc/source/whatsnew/index.rst @@ -21,7 +21,7 @@ Version 0.1.9 (31/10/2025) - :func:`~moocore.is_nondominated` and :func:`~moocore.filter_dominated` are faster for dimensions larger than 3. - ``moocore`` wheels are now built for ``aarch64`` (ARM64) in Linux and Windows. See the `installation instructions `_. -- :func:`~moocore.is_nondominated` and :func:`moocore.filter_dominated` are now stable in 2D and 3D with ``!keep_weakly``, that is, only the first of duplicated points is marked as nondominated. +- :func:`~moocore.is_nondominated` and :func:`~moocore.filter_dominated` are now stable in 2D and 3D with ``!keep_weakly``, that is, only the first of duplicated points is marked as nondominated. Version 0.1.8 (15/07/2025) From f40e169c1fab1b0fd77a5e03a643726ee3b1bd26 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 5 Nov 2025 19:41:35 +0000 Subject: [PATCH 19/49] * c/whv.c (print_point): Rename to debug_print_point. (print_rect): Rename to debug_print_rect. --- c/whv.c | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/c/whv.c b/c/whv.c index 029d9a0c7..7859ed3f6 100644 --- a/c/whv.c +++ b/c/whv.c @@ -75,12 +75,12 @@ double rect_weighted_hv2d(double *data, int n, double * rectangles, int rectangles_nrow, const double * reference) { -#define print_point(k, p, r, rect) \ - DEBUG2_PRINT("%d: p[%lu] = (%16.15g, %16.15g)" \ +#define debug_print_point(k, p, r, rect) \ + DEBUG2_PRINT("%d: p[%lu] = (%16.15g, %16.15g)" \ "\trectangle[%lu] = (%16.15g, %16.15g, %16.15g, %16.15g)\n", \ __LINE__, (unsigned long) k, p[0], p[1], (unsigned long) r, rect[0], rect[1], rect[2], rect[3]) -#define print_rect(r, rect) \ +#define debug_print_rect(r, rect) \ DEBUG2_PRINT("%d: rectangle[%lu] = (%16.15g, %16.15g, %16.15g, %16.15g, %16.15g)\n", \ __LINE__, (unsigned long) r, rect[0], rect[1], rect[2], rect[3], rect[4]) @@ -90,10 +90,10 @@ rect_weighted_hv2d(double *data, int n, double * rectangles, rect = rectangles + (ROW) * (nobj * 2 + 1); \ lower0= rect[0]; lower1= rect[1]; upper0= rect[2]; upper1= rect[3]; \ color = rect[4]; \ - print_rect(ROW, rect); \ - assert(lower0 < upper0); \ - assert(lower1 < upper1); \ - assert(color >= 0); \ + debug_print_rect(ROW, rect); \ + assert(lower0 < upper0); \ + assert(lower1 < upper1); \ + assert(color >= 0); \ } while(0) #define next_point() do { \ @@ -102,8 +102,9 @@ rect_weighted_hv2d(double *data, int n, double * rectangles, if (pk >= n || top == last_top || p[0] >= last_right) \ goto return_whv; \ p += nobj; \ - print_point(pk, p, r, rect); \ + debug_print_point(pk, p, r, rect); \ } while(0) + // We cannot use %zu for size_t because of MingW compiler. DEBUG2_PRINT("n = %lu\trectangles = %lu\n", (unsigned long)n, (unsigned long)rectangles_nrow); if (rectangles_nrow <= 0 || n <= 0) return 0; @@ -128,7 +129,7 @@ rect_weighted_hv2d(double *data, int n, double * rectangles, const double *p = data; int pk = 0; - print_point(pk, p, r, rect); + debug_print_point(pk, p, r, rect); double top = upper1; // lowest_upper1; const double last_top = rectangles[rectangles_nrow * (nobj * 2 + 1) - 2]; From cdde0c8b0f297bf8fa1d5d7d8b481b6f33d48d25 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 7 Nov 2025 15:52:02 +0000 Subject: [PATCH 20/49] c/hv_priv.h: Use restrict for x. --- c/hv_priv.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/c/hv_priv.h b/c/hv_priv.h index aa6569314..34c7e0b03 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -17,7 +17,7 @@ */ typedef struct dlnode { - const double * x; // point coordinates (objective vector). + const double * restrict x; // point coordinates (objective vector). struct dlnode * next[HV_DIMENSION - 2]; /* keeps the points sorted according to coordinates 2,3 and 4 (in the case of 2 and 3, only the points swept by 4 are kept) */ struct dlnode * prev[HV_DIMENSION - 2]; //keeps the points sorted according to coordinates 2 and 3 (except the sentinel 3) From 3654bd69d5612814239ff77f9910e7bc8cd61442 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 11 Nov 2025 12:03:27 +0000 Subject: [PATCH 21/49] c/hv4d.c: Remove unlikely. --- c/hv4d.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/c/hv4d.c b/c/hv4d.c index 87e0fe3ae..f9db3792d 100644 --- a/c/hv4d.c +++ b/c/hv4d.c @@ -132,7 +132,7 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n bool p_leq_new_2 = p_lt_new_2 | p_eq_new_2; // if (weakly_dominates(px, newx, 3)) { // Slower - if (unlikely(p_leq_new_0 & p_leq_new_1 & p_leq_new_2)) { + if (p_leq_new_0 & p_leq_new_1 & p_leq_new_2) { //new->ndomr++; assert(weakly_dominates(px, newx, 4)); return false; @@ -243,6 +243,8 @@ hv4dplusU(dlnode_t * list) add_to_z(new); update_links(list, new); } + // FIXME: It new was dominated, can we accumulate the height and update + // hv later? double height = new->next[1]->x[3] - new->x[3]; assert(height >= 0); hv += volume * height; From d570c484552dd9110ff442e27f68f93738e91dfe Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 12 Nov 2025 14:47:55 +0000 Subject: [PATCH 22/49] * c/gcc.mk: Add notes about GCC options. --- c/gcc.mk | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/c/gcc.mk b/c/gcc.mk index 1f11496f6..e68f3da09 100644 --- a/c/gcc.mk +++ b/c/gcc.mk @@ -7,7 +7,11 @@ WARN_CFLAGS = -pedantic -Wall -Wextra -Wvla -Wconversion -Wno-sign-conversion -W ifeq ($(DEBUG), 0) SANITIZERS ?= OPT_CFLAGS ?= -DNDEBUG -O3 -flto - # -ffinite-math-only -fno-signed-zeros allows further vectorization. + # GCC options that improve performance: + # -ffinite-math-only -fno-signed-zeros + # GCC options that may or may not improve performance: + # -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 OPT_CFLAGS ?= -g3 -O0 From 8b8e97d003c3229e38aa1a5845960da7b400ae32 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 5 Nov 2025 19:48:17 +0000 Subject: [PATCH 23/49] c/hv.c: Remove unused code. --- c/hv.c | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/c/hv.c b/c/hv.c index 10b6f8c7b..555d43e66 100644 --- a/c/hv.c +++ b/c/hv.c @@ -260,7 +260,8 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, /* Delete all points x[dim] > bound[d_stop]. In case of repeated coordinates, delete also all points x[dim] == bound[d_stop] except one. */ - while (p1->x[dim] > bound[d_stop] || p1->prev[d_stop]->x[dim] >= bound[d_stop]) { + while (p1->x[dim] > bound[d_stop] + || p1->prev[d_stop]->x[dim] >= bound[d_stop]) { // FIXME: Instead of deleting each point, unlink the start and end // nodes after the loop. delete(p1, dim, bound); @@ -279,11 +280,13 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, ASSUME(c == 1); update_area(p1->area, p1->x, ref, dim); p1->vol[d_stop] = 0; - if (p0->x == NULL) { + assert(p0->x != NULL); + /* if (p0->x == NULL) return p1->area[d_stop] * (ref[dim] - p1->x[dim]); - } + */ hyperv = p1->area[d_stop] * (p0->x[dim] - p1->x[dim]); - bound[d_stop] = p0->x[dim]; + // FIXME: This is never used. + // bound[d_stop] = p0->x[dim]; reinsert(p0, dim, bound); c++; p1 = p0; @@ -310,11 +313,13 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, } p1->area[d_stop] = hypera; if (p0->x == NULL) { + bound[d_stop] = p1->x[dim]; hyperv += hypera * (ref[dim] - p1->x[dim]); return hyperv; } hyperv += hypera * (p0->x[dim] - p1->x[dim]); - bound[d_stop] = p0->x[dim]; + // FIXME: This is never used. + // bound[d_stop] = p0->x[dim]; reinsert(p0, dim, bound); c++; p1 = p0; From 52d521b5e4502bc1023ae7c3aec4854d1823ec14 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 11 Nov 2025 12:20:22 +0000 Subject: [PATCH 24/49] c/hv.c (fpli_setup_cdllist): Simplify. --- c/hv.c | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/c/hv.c b/c/hv.c index 555d43e66..c99481cd2 100644 --- a/c/hv.c +++ b/c/hv.c @@ -73,8 +73,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, head->prev = head->next + d_stop * (n+1); head->area = malloc(2 * d_stop * (n+1) * sizeof(*data)); head->vol = head->area + d_stop * (n+1); - head->x = NULL; /* head contains no data */ - head->ignore = 0; /* should never get used */ + head->x = NULL; // head contains no data + head->ignore = 0; // should never get used size_t i = 1; for (size_t j = 0; j < n; j++) { @@ -82,7 +82,7 @@ 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. */ if (likely(strongly_dominates(data + j * d, ref, d))) { - head[i].x = data + (j+1) * d; /* this will be fixed a few lines below... */ + head[i].x = data + (j+1) * d; // this will be fixed a few lines below... head[i].ignore = 0; head[i].next = head->next + i * d_stop; head[i].prev = head->prev + i * d_stop; @@ -99,12 +99,11 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, for (i = 0; i < n; i++) scratch[i] = head + i + 1; - for (int k = d-1; k >= 0; k--) { + for (int j = d_stop - 1; j >= 0; j--) { + // We shift x because qsort cannot take the dimension to sort as an argument. for (i = 0; i < n; i++) scratch[i]->x--; - int j = k - STOP_DIMENSION; - if (j < 0) - continue; + // Sort each dimension independently. qsort(scratch, n, sizeof(*scratch), compare_node); head->next[j] = scratch[0]; scratch[0]->prev[j] = head; @@ -115,6 +114,9 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, scratch[n-1]->next[j] = head; head->prev[j] = scratch[n-1]; } + // Reset x to point to the first objective. + for (i = 0; i < n; i++) + scratch[i]->x -= STOP_DIMENSION; free(scratch); From da0f2f7ed91526a11bd3465c2cb407c4e02dfe86 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 6 Nov 2025 09:57:52 +0000 Subject: [PATCH 25/49] c/hv.c: Separate fpli_hv_ge5d() from hv_recursive(). Add debug counters. --- c/hv.c | 207 +++++++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 158 insertions(+), 49 deletions(-) diff --git a/c/hv.c b/c/hv.c index c99481cd2..05c780275 100644 --- a/c/hv.c +++ b/c/hv.c @@ -150,24 +150,37 @@ update_bound(double * restrict bound, const double * restrict x, dimension_t dim } static void -delete(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +delete_dom(fpli_dlnode_t * restrict nodep, dimension_t dim) { ASSUME(dim > STOP_DIMENSION); for (dimension_t d = 0; d < dim - STOP_DIMENSION; d++) { nodep->prev[d]->next[d] = nodep->next[d]; nodep->next[d]->prev[d] = nodep->prev[d]; } +} + +static void +delete(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +{ + delete_dom(nodep, dim); update_bound(bound, nodep->x, dim); } + static void -reinsert(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +reinsert_nobound(fpli_dlnode_t * restrict nodep, dimension_t dim) { ASSUME(dim > STOP_DIMENSION); for (dimension_t d = 0; d < dim - STOP_DIMENSION; d++) { nodep->prev[d]->next[d] = nodep; nodep->next[d]->prev[d] = nodep; } +} + +static void +reinsert(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +{ + reinsert_nobound(nodep, dim); update_bound(bound, nodep->x, dim); } @@ -207,15 +220,6 @@ fpli_hv4d_setup_cdllist(const fpli_dlnode_t * restrict pp, (list+2)->prev[d] = q; } -static double -one_point_hv(const double * restrict x, const double * restrict ref, dimension_t d) -{ - double hv = ref[0] - x[0]; - for (dimension_t i = 1; i < d; i++) - hv *= (ref[i] - x[i]); - return hv; -} - double hv4dplusU(dlnode_t * list); static double @@ -227,6 +231,15 @@ fpli_hv4d(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, size_t c) return hv; } +static double +one_point_hv(const double * restrict x, const double * restrict ref, dimension_t d) +{ + double hv = ref[0] - x[0]; + for (dimension_t i = 1; i < d; i++) + hv *= (ref[i] - x[i]); + return hv; +} + static inline void update_area(double * restrict area, const double * restrict x, const double * restrict ref, dimension_t dim) @@ -242,44 +255,54 @@ update_area(double * restrict area, const double * restrict x, area[d + 1] *= area[d]; } +#if DEBUG >= 1 +//#define HV_COUNTERS +static size_t debug_counter[6] = { 0 }; +#endif + static double hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, dimension_t dim, size_t c, const double * restrict ref, double * restrict bound) { - ASSUME(dim > STOP_DIMENSION); ASSUME(c > 1); + ASSUME(dim >= STOP_DIMENSION); + if (dim == STOP_DIMENSION) { + /*--------------------------------------- + base case of dimension 4 + --------------------------------------*/ + return fpli_hv4d(list, list4d, c); + } + ASSUME(dim > STOP_DIMENSION); /* ------------------------------------------------------ General case for dimensions higher than 4D ------------------------------------------------------ */ const dimension_t d_stop = dim - STOP_DIMENSION; - fpli_dlnode_t *p1 = list->prev[d_stop]; - for (fpli_dlnode_t *pp = p1; pp->x; pp = pp->prev[d_stop]) { + fpli_dlnode_t * p1 = list->prev[d_stop]; + for (fpli_dlnode_t * pp = p1; pp->x; pp = pp->prev[d_stop]) { if (pp->ignore < dim) pp->ignore = 0; } - fpli_dlnode_t *p0 = list; + fpli_dlnode_t * p1_prev = p1->prev[d_stop]; + fpli_dlnode_t * p0 = list; /* Delete all points x[dim] > bound[d_stop]. In case of repeated coordinates, delete also all points x[dim] == bound[d_stop] except one. */ while (p1->x[dim] > bound[d_stop] - || p1->prev[d_stop]->x[dim] >= bound[d_stop]) { + || p1_prev->x[dim] >= bound[d_stop]) { // FIXME: Instead of deleting each point, unlink the start and end // nodes after the loop. delete(p1, dim, bound); p0 = p1; p1 = p1->prev[d_stop]; + p1_prev = p1->prev[d_stop]; c--; if (c == 1) break; } - double hyperv = 0; - if (c > 1) { - hyperv = p1->prev[d_stop]->vol[d_stop] + p1->prev[d_stop]->area[d_stop] - * (p1->x[dim] - p1->prev[d_stop]->x[dim]); - } else { - ASSUME(c == 1); + double hyperv; + if (c == 1) { update_area(p1->area, p1->x, ref, dim); p1->vol[d_stop] = 0; assert(p0->x != NULL); @@ -287,31 +310,45 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, return p1->area[d_stop] * (ref[dim] - p1->x[dim]); */ hyperv = p1->area[d_stop] * (p0->x[dim] - p1->x[dim]); - // FIXME: This is never used. + // FIXME: This is never used? // bound[d_stop] = p0->x[dim]; reinsert(p0, dim, bound); c++; p1 = p0; + p1_prev = p0->prev[d_stop]; p0 = p0->next[d_stop]; + } else { + ASSUME(c > 1); + DEBUG1(debug_counter[0]++); + hyperv = p1_prev->vol[d_stop] + p1_prev->area[d_stop] + * (p1->x[dim] - p1_prev->x[dim]); + assert(p0 != p1_prev); + assert(p0 == p1->next[d_stop]); + // p0->x may be NULL here and thus we may return below. } assert(c > 1); while (true) { + // FIXME: This is not true in the first iteration if c > 1 previously. + // assert(p0 == p1->prev[d_stop]); + assert(p1_prev == p1->prev[d_stop]); p1->vol[d_stop] = hyperv; double hypera; if (p1->ignore >= dim) { - hypera = p1->prev[d_stop]->area[d_stop]; + DEBUG1(debug_counter[1]++); + hypera = p1_prev->area[d_stop]; } else { - if (dim - 1 == STOP_DIMENSION) { - /*--------------------------------------- - base case of dimension 4 - --------------------------------------*/ - hypera = fpli_hv4d(list, list4d, c); - } else { - hypera = hv_recursive(list, list4d, dim - 1, c, ref, bound); - } - if (hypera <= p1->prev[d_stop]->area[d_stop]) + hypera = hv_recursive(list, list4d, 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. */ + if (p1->ignore == dim - 1) { + DEBUG1(debug_counter[2]++); p1->ignore = dim; + } else if (hypera <= p1_prev->area[d_stop]) { + DEBUG1(debug_counter[3]++); + p1->ignore = dim; + } } p1->area[d_stop] = hypera; if (p0->x == NULL) { @@ -320,15 +357,96 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, return hyperv; } hyperv += hypera * (p0->x[dim] - p1->x[dim]); - // FIXME: This is never used. + // FIXME: This is never used? // bound[d_stop] = p0->x[dim]; reinsert(p0, dim, bound); c++; p1 = p0; + p1_prev = p0->prev[d_stop]; p0 = p0->next[d_stop]; } } +static double +fpli_hv_ge5d(fpli_dlnode_t * restrict list, dimension_t dim, size_t c, + const double * restrict ref) +{ + ASSUME(c > 1); + ASSUME(dim > STOP_DIMENSION); + const dimension_t d_stop = dim - STOP_DIMENSION; + ASSUME(0 < d_stop && d_stop < 255); // Silence -Walloc-size-larger-than= warning + double * bound = malloc(d_stop * sizeof(*bound)); + for (dimension_t i = 0; i < d_stop; i++) + bound[i] = -DBL_MAX; + dlnode_t * list4d = new_cdllist(c, ref); + + /* ------------------------------------------------------ + General case for dimensions higher than 4D + ------------------------------------------------------ */ + fpli_dlnode_t * p1 = list->prev[d_stop]; + // FIXME: This should be the initial state of the list when building it. + // Delete all points in dimensions < dim. + do { + delete_dom(p1, dim); + p1 = p1->prev[d_stop]; + c--; + } while (c > 1); + + update_area(p1->area, p1->x, ref, dim); + p1->vol[d_stop] = 0; + fpli_dlnode_t * p0 = p1->next[d_stop]; + assert(p0->x != NULL); + double hyperv = p1->area[d_stop] * (p0->x[dim] - p1->x[dim]); + // FIXME: This is never used? + // bound[d_stop] = p0->x[dim]; + reinsert_nobound(p0, dim); + p1 = p0; + fpli_dlnode_t * p1_prev = p0->prev[d_stop]; + p0 = p0->next[d_stop]; + c++; + + assert(c > 1); + while (true) { + // FIXME: This is not true in the first iteration if c > 1 previously. + //assert(p0 == p1->prev[d_stop]); + assert(p1_prev == p1->prev[d_stop]); + p1->vol[d_stop] = hyperv; + assert(p1->ignore == 0); + double hypera = hv_recursive(list, list4d, 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. */ + if (p1->ignore == dim - 1) { + DEBUG1(debug_counter[4]++); + p1->ignore = dim; + } else if (hypera <= p1_prev->area[d_stop]) { + DEBUG1(debug_counter[5]++); + p1->ignore = dim; + } + p1->area[d_stop] = hypera; + if (p0->x == NULL) { + free_cdllist(list4d); + free(bound); +#if defined(HV_COUNTERS) && DEBUG >= 1 + for (size_t i = 0; i < sizeof(debug_counter)/sizeof(size_t); i++) + fprintf(stderr, "debug_counter[%zu] = %zu\n", i, debug_counter[i]); +#endif + hyperv += hypera * (ref[dim] - p1->x[dim]); + return hyperv; + } + hyperv += hypera * (p0->x[dim] - p1->x[dim]); + // FIXME: This is never used? + // bound[d_stop] = p0->x[dim]; + // FIXME: Does updating the bound here matters? + reinsert(p0, dim, bound); + p1 = p0; + p1_prev = p0->prev[d_stop]; + p0 = p0->next[d_stop]; + c++; + } +} + + static double hv2d(const double * restrict data, size_t n, const double * restrict ref) { @@ -365,31 +483,22 @@ double fpli_hv(const double * restrict data, int d, int npoints, { size_t n = (size_t) npoints; if (unlikely(n == 0)) return 0.0; - ASSUME(d < 256); - ASSUME(d > 1); + ASSUME(d > 1 && d < 256); if (d == 4) return hv4d(data, n, ref); if (d == 3) return hv3d(data, n, ref); if (d == 2) return hv2d(data, n, ref); dimension_t dim = (dimension_t) d; fpli_dlnode_t * list = fpli_setup_cdllist(data, dim, &n, ref); double hyperv; - if (unlikely(n == 0)) { - /* Returning here would leak memory. */ - hyperv = 0.0; + if (likely(n > 1)) { + hyperv = fpli_hv_ge5d(list, dim - 1, n, ref); } else if (unlikely(n == 1)) { hyperv = one_point_hv(list->next[0]->x, ref, dim); } else { - const dimension_t d_stop = dim - STOP_DIMENSION; - ASSUME(d_stop > 1 && d_stop < 255); // Silence -Walloc-size-larger-than= warning - double * bound = malloc(d_stop * sizeof(*bound)); - for (dimension_t i = 0; i < d_stop; i++) - bound[i] = -DBL_MAX; - dlnode_t * list4d = new_cdllist(n, ref); - hyperv = hv_recursive(list, list4d, dim - 1, n, ref, bound); - free_cdllist(list4d); - free(bound); + assert(n == 0); + hyperv = 0.0; // Returning here would leak memory. } - /* Clean up. */ + // Clean up. fpli_free_cdllist(list); return hyperv; } From ccba3207b85a1edd31d1d7ba624a51eedc287e8a Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 13 Nov 2025 16:23:13 +0000 Subject: [PATCH 26/49] * c/gcc_attribs.h (ASAN_POISON_MEMORY_REGION): New. --- c/gcc_attribs.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/c/gcc_attribs.h b/c/gcc_attribs.h index f69ab12c2..00418aa34 100644 --- a/c/gcc_attribs.h +++ b/c/gcc_attribs.h @@ -227,4 +227,10 @@ #define _attr_optimize_finite_math /* nothing */ #endif +#ifdef __SANITIZE_ADDRESS__ +# include +#else +# define ASAN_POISON_MEMORY_REGION(addr, size) ((void) (addr), (void) (size)) +#endif + #endif /* GCC_ATTRIBUTES */ From 7eef75e136510448f033baec4c1d0f159c0ed3dc Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 13 Nov 2025 15:41:27 +0000 Subject: [PATCH 27/49] * c/hv.c: Merge fpli_dlnode_t into dlnode_t. --- c/Makefile | 1 + c/hv.c | 242 +++++++++++++++++++++++------------------------ c/hv4d.c | 255 +------------------------------------------------- c/hv4d_priv.h | 248 ++++++++++++++++++++++++++++++++++++++++++++++++ c/hv_priv.h | 64 ++++++++----- c/libhv.mk | 2 +- 6 files changed, 411 insertions(+), 401 deletions(-) create mode 100644 c/hv4d_priv.h diff --git a/c/Makefile b/c/Makefile index e8dab4a06..e72cffba4 100644 --- a/c/Makefile +++ b/c/Makefile @@ -96,6 +96,7 @@ HDRS = avl.h \ hvapprox.h \ hv_priv.h \ hv3d_priv.h \ + hv4d_priv.h \ igd.h \ io.h \ io_priv.h \ diff --git a/c/hv.c b/c/hv.c index 05c780275..0e4abac62 100644 --- a/c/hv.c +++ b/c/hv.c @@ -34,48 +34,53 @@ #include #include "common.h" #include "hv.h" -#define HV_DIMENSION 4 -#include "hv_priv.h" +#define HV_RECURSIVE +#include "hv4d_priv.h" #define STOP_DIMENSION 3 // default: stop on dimension 4. -typedef struct fpli_dlnode { - const double * restrict x; // point vector - struct fpli_dlnode ** next; // next-node vector - struct fpli_dlnode ** prev; // previous-node vector - double * restrict area; // partial area - double * restrict vol; // partial volume - dimension_t ignore; // [0, 255] -} fpli_dlnode_t; - - static int compare_node(const void * restrict p1, const void * restrict p2) { - const double * restrict x1 = (*(const fpli_dlnode_t **)p1)->x; - const double * restrict x2 = (*(const fpli_dlnode_t **)p2)->x; + const double * restrict x1 = (*(const dlnode_t **)p1)->x; + const double * restrict x2 = (*(const dlnode_t **)p2)->x; return cmp_double_asc(*x1, *x2); } -/* - * Setup circular double-linked list in each dimension - */ +/** Setup circular double-linked list in each dimension. + + There are in fact two separate lists that are keep in sync: + + - A multi-dimensional list for dimensions 5 and above tracked by ->r_next[] and ->r_prev[], where head is the sentinel with head->x == NULL. -static fpli_dlnode_t * + - A list for dimensions 3 and 4 tracked by ->next[0 or 1] ->prev[0 or 1]. This list has 3 sentinels as required by hv4dplusU(). The first sentinel is saved in head->next[0]. + + */ +static dlnode_t * fpli_setup_cdllist(const double * restrict data, dimension_t d, size_t * restrict size, const double * restrict ref) { - ASSUME(d > STOP_DIMENSION); + ASSUME(d > STOP_DIMENSION + 1); dimension_t d_stop = d - STOP_DIMENSION; size_t n = *size; - fpli_dlnode_t * head = malloc((n+1) * sizeof(*head)); + + dlnode_t * head = malloc((n+1) * sizeof(*head)); // Allocate single blocks of memory as much as possible. - head->next = malloc(2 * d_stop * (n+1) * sizeof(head)); - head->prev = head->next + d_stop * (n+1); + // We need space in r_next and r_prev for dimension 5 and above (d_stop - 1). + head->r_next = malloc(2 * (d_stop - 1) * (n+1) * sizeof(head)); + head->r_prev = head->r_next + (d_stop - 1) * (n+1); + // We only need space in area and vol for dimension 4 and above. head->area = malloc(2 * d_stop * (n+1) * sizeof(*data)); head->vol = head->area + d_stop * (n+1); head->x = NULL; // head contains no data head->ignore = 0; // should never get used + // Reserve space for the sentinels. + dlnode_t * list4d = new_cdllist(0, ref); + // 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. + size_t i = 1; for (size_t j = 0; j < n; j++) { /* Filters those points that do not strictly dominate the reference @@ -84,8 +89,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, if (likely(strongly_dominates(data + j * d, ref, d))) { head[i].x = data + (j+1) * d; // this will be fixed a few lines below... head[i].ignore = 0; - head[i].next = head->next + i * d_stop; - head[i].prev = head->prev + i * d_stop; + 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; i++; @@ -95,24 +100,39 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, if (unlikely(n == 0)) goto finish; - fpli_dlnode_t **scratch = malloc(n * sizeof(*scratch)); + dlnode_t **scratch = malloc(n * sizeof(*scratch)); for (i = 0; i < n; i++) scratch[i] = head + i + 1; - for (int j = d_stop - 1; j >= 0; j--) { - // We shift x because qsort cannot take the dimension to sort as an argument. + for (int j = d_stop - 2; j >= -1; j--) { + /* FIXME: replace qsort() by something better: + https://github.com/numpy/x86-simd-sort + https://github.com/google/highway/tree/52a2d98d07852c5d69284e175666e5f8cc7d8285/hwy/contrib/sort + */ + // We shift x because qsort() cannot take the dimension to sort as an argument. for (i = 0; i < n; i++) scratch[i]->x--; // Sort each dimension independently. qsort(scratch, n, sizeof(*scratch), compare_node); - head->next[j] = scratch[0]; - scratch[0]->prev[j] = head; - for (i = 1; i < n; i++) { - scratch[i-1]->next[j] = scratch[i]; - scratch[i]->prev[j] = scratch[i-1]; + if (j == -1) { + (list4d+1)->next[1] = scratch[0]; + scratch[0]->prev[1] = list4d+1; + for (i = 1; i < n; i++) { + scratch[i-1]->next[1] = scratch[i]; + scratch[i]->prev[1] = scratch[i-1]; + } + scratch[n-1]->next[1] = list4d+2; + (list4d+2)->prev[1] = scratch[n-1]; + } 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]->next[j] = head; - head->prev[j] = scratch[n-1]; } // Reset x to point to the first objective. for (i = 0; i < n; i++) @@ -120,18 +140,20 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, free(scratch); - // FIXME: This should not be necessary. - for (i = 0; i < d_stop; i++) - head->area[i] = 0; + // 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; } -static void fpli_free_cdllist(fpli_dlnode_t * head) +static void fpli_free_cdllist(dlnode_t * head) { - free(head->next); + assert(head->next[0] == head->prev[0]); + free_cdllist(head->next[0]); // free 4D sentinels + free(head->r_next); free(head->area); free(head); } @@ -150,17 +172,22 @@ update_bound(double * restrict bound, const double * restrict x, dimension_t dim } static void -delete_dom(fpli_dlnode_t * restrict nodep, dimension_t dim) +delete_dom(dlnode_t * restrict nodep, dimension_t dim) { ASSUME(dim > STOP_DIMENSION); - for (dimension_t d = 0; d < dim - STOP_DIMENSION; d++) { - nodep->prev[d]->next[d] = nodep->next[d]; - nodep->next[d]->prev[d] = nodep->prev[d]; + assert(nodep->x); + // d=0 is dimension 5. + for (dimension_t d = 0; d < dim - 1 - STOP_DIMENSION; d++) { + 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]; } static void -delete(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +delete(dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) { delete_dom(nodep, dim); update_bound(bound, nodep->x, dim); @@ -168,65 +195,35 @@ delete(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) static void -reinsert_nobound(fpli_dlnode_t * restrict nodep, dimension_t dim) +reinsert_nobound(dlnode_t * restrict nodep, dimension_t dim) { ASSUME(dim > STOP_DIMENSION); - for (dimension_t d = 0; d < dim - STOP_DIMENSION; d++) { - nodep->prev[d]->next[d] = nodep; - nodep->next[d]->prev[d] = nodep; + assert(nodep->x); + // d=0 is dimension 5. + for (dimension_t d = 0; d < dim - 1 - STOP_DIMENSION; d++) { + 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; } static void -reinsert(fpli_dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) +reinsert(dlnode_t * restrict nodep, dimension_t dim, double * restrict bound) { reinsert_nobound(nodep, dim); update_bound(bound, nodep->x, dim); } -static void -fpli_hv4d_setup_cdllist(const fpli_dlnode_t * restrict pp, - dlnode_t * restrict list, size_t n _attr_maybe_unused) -{ - ASSUME(n > 1); - reset_sentinels(list); - - const dimension_t d = HV_DIMENSION - 3; // index within the list. - dlnode_t * q = list+1; - dlnode_t * list3 = list+3; - assert(list->next[d] == list + 1); - assert(q->next[d] == list + 2); - for (size_t i = 0; pp->x != NULL; pp = pp->next[0]) { - dlnode_t * p = list3 + i; - p->x = pp->x; - // Initialize it when debugging so it will crash if uninitialized. - DEBUG1( - p->closest[0] = NULL; - p->closest[1] = NULL; - p->cnext[0] = NULL; - p->cnext[1] = NULL;); - // FIXME: Can we use pp->ignore to initialize p->ndomr? - //p->ndomr = 0; - // Link the list in order. - q->next[d] = p; - p->prev[d] = q; - q = p; - i++; - } - assert((list3 + n - 1) == q); - assert(list+2 == list->prev[d]); - // q = last point, q->next = s3, s3->prev = last point - q->next[d] = list+2; - (list+2)->prev[d] = q; -} - -double hv4dplusU(dlnode_t * list); - static double -fpli_hv4d(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, size_t c) +fpli_hv4d(dlnode_t * restrict list, size_t c _attr_maybe_unused) { ASSUME(c > 1); - fpli_hv4d_setup_cdllist(list->next[0], list4d, c); + assert(list->next[0] == list->prev[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; } @@ -255,14 +252,11 @@ update_area(double * restrict area, const double * restrict x, area[d + 1] *= area[d]; } -#if DEBUG >= 1 //#define HV_COUNTERS -static size_t debug_counter[6] = { 0 }; -#endif +_attr_maybe_unused static size_t debug_counter[6] = { 0 }; static double -hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, - dimension_t dim, size_t c, +hv_recursive(dlnode_t * restrict list, dimension_t dim, size_t c, const double * restrict ref, double * restrict bound) { ASSUME(c > 1); @@ -271,20 +265,22 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, /*--------------------------------------- base case of dimension 4 --------------------------------------*/ - return fpli_hv4d(list, list4d, c); + return fpli_hv4d(list, c); } ASSUME(dim > STOP_DIMENSION); /* ------------------------------------------------------ General case for dimensions higher than 4D ------------------------------------------------------ */ const dimension_t d_stop = dim - STOP_DIMENSION; - fpli_dlnode_t * p1 = list->prev[d_stop]; - for (fpli_dlnode_t * pp = p1; pp->x; pp = pp->prev[d_stop]) { + assert(d_stop > 0); + // d_stop - 1 is dimension 5 in r_prev and r_next. + dlnode_t * p1 = list->r_prev[d_stop - 1]; + for (dlnode_t * pp = p1; pp->x; pp = pp->r_prev[d_stop - 1]) { if (pp->ignore < dim) pp->ignore = 0; } - fpli_dlnode_t * p1_prev = p1->prev[d_stop]; - fpli_dlnode_t * p0 = list; + dlnode_t * p1_prev = p1->r_prev[d_stop - 1]; + dlnode_t * p0 = list; /* Delete all points x[dim] > bound[d_stop]. In case of repeated coordinates, delete also all points x[dim] == bound[d_stop] except one. */ @@ -294,8 +290,8 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, // nodes after the loop. delete(p1, dim, bound); p0 = p1; - p1 = p1->prev[d_stop]; - p1_prev = p1->prev[d_stop]; + p1 = p1->r_prev[d_stop - 1]; + p1_prev = p1->r_prev[d_stop - 1]; c--; if (c == 1) break; @@ -315,30 +311,30 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, reinsert(p0, dim, bound); c++; p1 = p0; - p1_prev = p0->prev[d_stop]; - p0 = p0->next[d_stop]; + p1_prev = p0->r_prev[d_stop - 1]; + p0 = p0->r_next[d_stop - 1]; } else { ASSUME(c > 1); DEBUG1(debug_counter[0]++); hyperv = p1_prev->vol[d_stop] + p1_prev->area[d_stop] * (p1->x[dim] - p1_prev->x[dim]); assert(p0 != p1_prev); - assert(p0 == p1->next[d_stop]); + assert(p0 == p1->r_next[d_stop - 1]); // p0->x may be NULL here and thus we may return below. } assert(c > 1); while (true) { // FIXME: This is not true in the first iteration if c > 1 previously. - // assert(p0 == p1->prev[d_stop]); - assert(p1_prev == p1->prev[d_stop]); + // assert(p0 == p1->r_prev[d_stop - 1]); + assert(p1_prev == p1->r_prev[d_stop - 1]); p1->vol[d_stop] = hyperv; double hypera; if (p1->ignore >= dim) { DEBUG1(debug_counter[1]++); hypera = p1_prev->area[d_stop]; } else { - hypera = hv_recursive(list, list4d, dim - 1, c, ref, bound); + 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. */ @@ -362,13 +358,13 @@ hv_recursive(fpli_dlnode_t * restrict list, dlnode_t * restrict list4d, reinsert(p0, dim, bound); c++; p1 = p0; - p1_prev = p0->prev[d_stop]; - p0 = p0->next[d_stop]; + p1_prev = p0->r_prev[d_stop - 1]; + p0 = p0->r_next[d_stop - 1]; } } static double -fpli_hv_ge5d(fpli_dlnode_t * restrict list, dimension_t dim, size_t c, +fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, const double * restrict ref) { ASSUME(c > 1); @@ -378,41 +374,40 @@ fpli_hv_ge5d(fpli_dlnode_t * restrict list, dimension_t dim, size_t c, double * bound = malloc(d_stop * sizeof(*bound)); for (dimension_t i = 0; i < d_stop; i++) bound[i] = -DBL_MAX; - dlnode_t * list4d = new_cdllist(c, ref); /* ------------------------------------------------------ General case for dimensions higher than 4D ------------------------------------------------------ */ - fpli_dlnode_t * p1 = list->prev[d_stop]; + dlnode_t * p1 = list->r_prev[d_stop - 1]; // FIXME: This should be the initial state of the list when building it. // Delete all points in dimensions < dim. do { delete_dom(p1, dim); - p1 = p1->prev[d_stop]; + p1 = p1->r_prev[d_stop - 1]; c--; } while (c > 1); update_area(p1->area, p1->x, ref, dim); p1->vol[d_stop] = 0; - fpli_dlnode_t * p0 = p1->next[d_stop]; + dlnode_t * p0 = p1->r_next[d_stop - 1]; assert(p0->x != NULL); double hyperv = p1->area[d_stop] * (p0->x[dim] - p1->x[dim]); // FIXME: This is never used? // bound[d_stop] = p0->x[dim]; reinsert_nobound(p0, dim); p1 = p0; - fpli_dlnode_t * p1_prev = p0->prev[d_stop]; - p0 = p0->next[d_stop]; + 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->prev[d_stop]); - assert(p1_prev == p1->prev[d_stop]); + //assert(p0 == p1->r_prev[d_stop]); + assert(p1_prev == p1->r_prev[d_stop - 1]); p1->vol[d_stop] = hyperv; assert(p1->ignore == 0); - double hypera = hv_recursive(list, list4d, dim - 1, c, ref, bound); + double 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. */ @@ -425,7 +420,6 @@ fpli_hv_ge5d(fpli_dlnode_t * restrict list, dimension_t dim, size_t c, } p1->area[d_stop] = hypera; if (p0->x == NULL) { - free_cdllist(list4d); free(bound); #if defined(HV_COUNTERS) && DEBUG >= 1 for (size_t i = 0; i < sizeof(debug_counter)/sizeof(size_t); i++) @@ -440,8 +434,8 @@ fpli_hv_ge5d(fpli_dlnode_t * restrict list, dimension_t dim, size_t c, // FIXME: Does updating the bound here matters? reinsert(p0, dim, bound); p1 = p0; - p1_prev = p0->prev[d_stop]; - p0 = p0->next[d_stop]; + p1_prev = p0->r_prev[d_stop - 1]; + p0 = p0->r_next[d_stop - 1]; c++; } } @@ -488,12 +482,12 @@ double fpli_hv(const double * restrict data, int d, int npoints, if (d == 3) return hv3d(data, n, ref); if (d == 2) return hv2d(data, n, ref); dimension_t dim = (dimension_t) d; - fpli_dlnode_t * list = fpli_setup_cdllist(data, dim, &n, ref); + dlnode_t * list = fpli_setup_cdllist(data, dim, &n, ref); double hyperv; if (likely(n > 1)) { hyperv = fpli_hv_ge5d(list, dim - 1, n, ref); } else if (unlikely(n == 1)) { - hyperv = one_point_hv(list->next[0]->x, ref, dim); + hyperv = one_point_hv(list->r_next[0]->x, ref, dim); } else { assert(n == 0); hyperv = 0.0; // Returning here would leak memory. diff --git a/c/hv4d.c b/c/hv4d.c index f9db3792d..570aa4af4 100644 --- a/c/hv4d.c +++ b/c/hv4d.c @@ -1,257 +1,4 @@ -/****************************************************************************** - HV4D+ algorithm. - ------------------------------------------------------------------------------ - - Copyright (C) 2013, 2016, 2017 - 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. - - ------------------------------------------------------------------------------ - - Modified by Manuel López-Ibáñez (2025): - - Integration within moocore. - - Correct handling of weakly dominated points and repeated coordinates during - preprocessing(). - - More efficient setup_cdllist() and preprocessing() in terms of time and memory. - -******************************************************************************/ - -#include -#include -#include "common.h" -#define HV_DIMENSION 4 -#include "hv_priv.h" - -// ---------- Data Structures Functions --------------------------------------- - -static inline void _attr_maybe_unused -print_point(const char *s, const double * x) -{ - fprintf(stderr, "%s: %g %g %g %g\n", s, x[0], x[1], x[2], x[4]); -} - - - - -// ------------ Update data structure ----------------------------------------- - -static inline void -add_to_z(dlnode_t * new) -{ - new->next[0] = new->prev[0]->next[0]; //in case new->next[0] was removed for being dominated - new->next[0]->prev[0] = new; - new->prev[0]->next[0] = new; -} - -static inline bool -lex_cmp_3d_102(const double * restrict a, const double * restrict b) -{ - return a[1] < b[1] || (a[1] == b[1] && (a[0] < b[0] || (a[0] == b[0] && a[2] < b[2]))); -} - -static inline bool -lex_cmp_3d_012(const double * restrict a, const double * restrict b) -{ - return a[0] < b[0] || (a[0] == b[0] && (a[1] < b[1] || (a[1] == b[1] && a[2] < b[2]))); -} - -/* - Go through the points in the order of z and either remove points that are - dominated by new with respect to x,y,z or update the cx and cy lists by - adding new. -*/ -static void -update_links(dlnode_t * restrict list, dlnode_t * restrict new) -{ - assert(list+2 == list->prev[0]); - const double newx[] = { new->x[0], new->x[1], new->x[2] }; - dlnode_t * p = new->next[0]; - const dlnode_t * stop = list+2; - while (p != stop) { - const double * px = p->x; - // px dominates newx (but not equal) - if (px[0] <= newx[0] && px[1] <= newx[1] && (px[0] < newx[0] || px[1] < newx[1])) - return; - - if (newx[0] <= px[0]){ - //new <= p - if (newx[1] <= px[1]){ - assert(weakly_dominates(newx, px, 3)); - //p->ndomr++; - remove_from_z(p); - } else if (newx[0] < px[0] && lex_cmp_3d_102(newx, p->closest[1]->x)) { // newx[1] > px[1] - p->closest[1] = new; - } - } else if (newx[1] < px[1] && lex_cmp_3d_012(newx, p->closest[0]->x)) { //newx[0] > px[0] - p->closest[0] = new; - } - p = p->next[0]; - } -} - - - -// This does what setupZandClosest does while reconstructing L at z = new->x[2]. -_attr_optimize_finite_math -__attribute__((hot)) -static bool -restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict new) -{ - // FIXME: This is the most expensive function in the HV4D+ algorithm. - const double newx[] = { new->x[0], new->x[1], new->x[2], new->x[3] }; - assert(list+1 == list->next[0]); - dlnode_t * closest0 = list+1; - dlnode_t * closest1 = list; - const double * closest0x = closest0->x; - const double * closest1x = closest1->x; - dlnode_t * p = (list+1)->next[0]; - assert(p == list->next[0]->next[0]); - restart_list_y(list); - while (true) { - const double * restrict px = p->x; - // Help auto-vectorization. - bool p_lt_new_0 = px[0] < newx[0]; - bool p_lt_new_1 = px[1] < newx[1]; - bool p_lt_new_2 = px[2] < newx[2]; - bool p_eq_new_0 = px[0] == newx[0]; - bool p_eq_new_1 = px[1] == newx[1]; - bool p_eq_new_2 = px[2] == newx[2]; - bool p_leq_new_0 = p_lt_new_0 | p_eq_new_0; - bool p_leq_new_1 = p_lt_new_1 | p_eq_new_1; - bool p_leq_new_2 = p_lt_new_2 | p_eq_new_2; - - // if (weakly_dominates(px, newx, 3)) { // Slower - if (p_leq_new_0 & p_leq_new_1 & p_leq_new_2) { - //new->ndomr++; - assert(weakly_dominates(px, newx, 4)); - 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))))) { - new->closest[0] = closest0; - new->closest[1] = closest1; - new->prev[0] = p->prev[0]; - new->next[0] = p; - return true; - } - - // FIXME: Can we move reconstruct() after setup_z_and_closest() ? - // reconstruct() - set_cnext_to_closest(p); - p->cnext[0]->cnext[1] = p; - p->cnext[1]->cnext[0] = p; - - // setup_z_and_closest() - ASSUME(!p_leq_new_0 || !p_leq_new_1); - if (p_lt_new_1 && (px[0] < closest0x[0] || (px[0] == closest0x[0] && px[1] < closest0x[1]))) { - closest0 = p; - closest0x = px; - } else if (p_lt_new_0 && (px[1] < closest1x[1] || (px[1] == closest1x[1] && px[0] < closest1x[0]))) { - closest1 = p; - closest1x = px; - } - p = p->next[0]; - } - unreachable(); -} - -// 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 new) -{ - set_cnext_to_closest(new); - const double * newx = new->x; - // if newx[0] == new->cnext[0]->x[0], the first area is zero - double area = compute_area_no_inners(newx, new->cnext[0], 1); - double volume = 0; - double lastz = newx[2]; - dlnode_t * p = new->next[0]; - assert(!weakly_dominates(p->x, newx, 4)); - while (true) { - const double * px = p->x; - volume += area * (px[2] - lastz); - - if (px[0] <= newx[0] && px[1] <= newx[1]) - return volume; - - ASSUME(px[0] > newx[0] || px[1] > newx[1]); - assert(!weakly_dominates(px, p->next[0]->x, 4)); - - set_cnext_to_closest(p); - - if (px[0] < newx[0]) { - if (px[1] <= new->cnext[1]->x[1]) { - const double tmpx[] = { newx[0], px[1] }; - // if px[1] == new->cnext[1]->x[1] then area starts at 0. - area -= compute_area_no_inners(tmpx, new->cnext[1], 0); - p->cnext[1] = new->cnext[1]; - p->cnext[0]->cnext[1] = p; - new->cnext[1] = p; - } - } else if (px[1] < newx[1]) { - if (px[0] <= new->cnext[0]->x[0]) { - const double tmpx[] = { px[0], newx[1] }; - // if px[0] == new->cnext[0]->x[0] then area starts at 0. - area -= compute_area_no_inners(tmpx, new->cnext[0], 1); - p->cnext[0] = new->cnext[0]; - p->cnext[1]->cnext[0] = p; - new->cnext[0] = p; - } - } else { - assert(px[0] >= newx[0] && px[1] >= newx[1]); - // if px[0] == p->cnext[0]->x[0] then area starts at 0. - area -= compute_area_no_inners(px, p->cnext[0], 1); - p->cnext[1]->cnext[0] = p; - p->cnext[0]->cnext[1] = p; - } - lastz = px[2]; - p = p->next[0]; - } - unreachable(); -} - -/* Compute the hypervolume indicator in d=4 by iteratively computing the one - contribution problem in d=3. */ -double -hv4dplusU(dlnode_t * list) -{ - assert(list+2 == list->prev[0]); - assert(list+2 == list->prev[1]); - assert(list+1 == list->next[1]); - - double volume = 0, hv = 0; - dlnode_t * new = (list+1)->next[1]; - dlnode_t * last = list+2; - while (new != last) { - if (restart_base_setup_z_and_closest(list, new)) { - // new was not dominated by something else. - double new_v = one_contribution_3d(new); - assert(new_v > 0); - volume += new_v; - add_to_z(new); - update_links(list, new); - } - // FIXME: It new was dominated, can we accumulate the height and update - // hv later? - double height = new->next[1]->x[3] - new->x[3]; - assert(height >= 0); - hv += volume * height; - new = new->next[1]; - } - return hv; -} +#include "hv4d_priv.h" double hv4d(const double * restrict data, size_t n, const double * restrict ref) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h new file mode 100644 index 000000000..0adaddfc0 --- /dev/null +++ b/c/hv4d_priv.h @@ -0,0 +1,248 @@ +#ifndef _HV4D_PRIV_H +#define _HV4D_PRIV_H + +/****************************************************************************** + HV4D+ algorithm. + ------------------------------------------------------------------------------ + + Copyright (C) 2013, 2016, 2017 + 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. + + ------------------------------------------------------------------------------ + + Modified by Manuel López-Ibáñez (2025): + - Integration within moocore. + - Correct handling of weakly dominated points and repeated coordinates during + preprocessing(). + - More efficient setup_cdllist() and preprocessing() in terms of time and memory. + +******************************************************************************/ + +#include +#include +#include "common.h" +#define HV_DIMENSION 4 +#include "hv_priv.h" + +// ------------ Update data structure ----------------------------------------- + +static inline void +add_to_z(dlnode_t * newp) +{ + newp->next[0] = newp->prev[0]->next[0]; //in case newp->next[0] was removed for being dominated + newp->next[0]->prev[0] = newp; + newp->prev[0]->next[0] = newp; +} + +static inline bool +lex_cmp_3d_102(const double * restrict a, const double * restrict b) +{ + return a[1] < b[1] || (a[1] == b[1] && (a[0] < b[0] || (a[0] == b[0] && a[2] < b[2]))); +} + +static inline bool +lex_cmp_3d_012(const double * restrict a, const double * restrict b) +{ + return a[0] < b[0] || (a[0] == b[0] && (a[1] < b[1] || (a[1] == b[1] && a[2] < b[2]))); +} + +/* + Go through the points in the order of z and either remove points that are + dominated by newp with respect to x,y,z or update the cx and cy lists by + adding newp. +*/ +static void +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; + while (p != stop) { + const double * px = p->x; + // px dominates newx (but not equal) + if (px[0] <= newx[0] && px[1] <= newx[1] && (px[0] < newx[0] || px[1] < newx[1])) + return; + + if (newx[0] <= px[0]){ + //new <= p + if (newx[1] <= px[1]){ + assert(weakly_dominates(newx, px, 3)); + //p->ndomr++; + remove_from_z(p); + } else if (newx[0] < px[0] && lex_cmp_3d_102(newx, p->closest[1]->x)) { // newx[1] > px[1] + p->closest[1] = newp; + } + } else if (newx[1] < px[1] && lex_cmp_3d_012(newx, p->closest[0]->x)) { //newx[0] > px[0] + p->closest[0] = newp; + } + p = p->next[0]; + } +} + + + +// This does what setupZandClosest does while reconstructing L at z = newp->x[2]. +_attr_optimize_finite_math +__attribute__((hot)) +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. + const double newx[] = { newp->x[0], newp->x[1], newp->x[2], newp->x[3] }; + assert(list+1 == list->next[0]); + dlnode_t * closest0 = list+1; + dlnode_t * closest1 = list; + const double * closest0x = closest0->x; + const double * closest1x = closest1->x; + dlnode_t * p = (list+1)->next[0]; + assert(p == list->next[0]->next[0]); + restart_list_y(list); + while (true) { + const double * restrict px = p->x; + // Help auto-vectorization. + bool p_lt_new_0 = px[0] < newx[0]; + bool p_lt_new_1 = px[1] < newx[1]; + bool p_lt_new_2 = px[2] < newx[2]; + bool p_eq_new_0 = px[0] == newx[0]; + bool p_eq_new_1 = px[1] == newx[1]; + bool p_eq_new_2 = px[2] == newx[2]; + bool p_leq_new_0 = p_lt_new_0 | p_eq_new_0; + bool p_leq_new_1 = p_lt_new_1 | p_eq_new_1; + bool p_leq_new_2 = p_lt_new_2 | p_eq_new_2; + + // if (weakly_dominates(px, newx, 3)) { // Slower + if (p_leq_new_0 & p_leq_new_1 & p_leq_new_2) { + //new->ndomr++; + assert(weakly_dominates(px, newx, 4)); + 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))))) { + newp->closest[0] = closest0; + newp->closest[1] = closest1; + newp->prev[0] = p->prev[0]; + newp->next[0] = p; + return true; + } + + // FIXME: Can we move reconstruct() after setup_z_and_closest() ? + // reconstruct() + set_cnext_to_closest(p); + p->cnext[0]->cnext[1] = p; + p->cnext[1]->cnext[0] = p; + + // setup_z_and_closest() + ASSUME(!p_leq_new_0 || !p_leq_new_1); + if (p_lt_new_1 && (px[0] < closest0x[0] || (px[0] == closest0x[0] && px[1] < closest0x[1]))) { + closest0 = p; + closest0x = px; + } else if (p_lt_new_0 && (px[1] < closest1x[1] || (px[1] == closest1x[1] && px[0] < closest1x[0]))) { + closest1 = p; + closest1x = px; + } + p = p->next[0]; + } + unreachable(); +} + +// 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) +{ + set_cnext_to_closest(newp); + const double * newx = newp->x; + // if newx[0] == newp->cnext[0]->x[0], the first area is zero + double area = compute_area_no_inners(newx, newp->cnext[0], 1); + double volume = 0; + double lastz = newx[2]; + dlnode_t * p = newp->next[0]; + assert(!weakly_dominates(p->x, newx, 4)); + while (true) { + const double * px = p->x; + volume += area * (px[2] - lastz); + + if (px[0] <= newx[0] && px[1] <= newx[1]) + return volume; + + ASSUME(px[0] > newx[0] || px[1] > newx[1]); + assert(!weakly_dominates(px, p->next[0]->x, 4)); + + set_cnext_to_closest(p); + + if (px[0] < newx[0]) { + if (px[1] <= newp->cnext[1]->x[1]) { + const double tmpx[] = { newx[0], px[1] }; + // if px[1] == newp->cnext[1]->x[1] then area starts at 0. + area -= compute_area_no_inners(tmpx, newp->cnext[1], 0); + p->cnext[1] = newp->cnext[1]; + p->cnext[0]->cnext[1] = p; + newp->cnext[1] = p; + } + } else if (px[1] < newx[1]) { + if (px[0] <= newp->cnext[0]->x[0]) { + const double tmpx[] = { px[0], newx[1] }; + // if px[0] == newp->cnext[0]->x[0] then area starts at 0. + area -= compute_area_no_inners(tmpx, newp->cnext[0], 1); + p->cnext[0] = newp->cnext[0]; + p->cnext[1]->cnext[0] = p; + newp->cnext[0] = p; + } + } else { + assert(px[0] >= newx[0] && px[1] >= newx[1]); + // if px[0] == p->cnext[0]->x[0] then area starts at 0. + area -= compute_area_no_inners(px, p->cnext[0], 1); + p->cnext[1]->cnext[0] = p; + p->cnext[0]->cnext[1] = p; + } + lastz = px[2]; + p = p->next[0]; + } + unreachable(); +} + +/* Compute the hypervolume indicator in d=4 by iteratively computing the one + contribution problem in d=3. */ +static double +hv4dplusU(dlnode_t * list) +{ + assert(list+2 == list->prev[0]); + assert(list+2 == list->prev[1]); + assert(list+1 == list->next[1]); + + double volume = 0, hv = 0; + dlnode_t * newp = (list+1)->next[1]; + const dlnode_t * last = list+2; + while (newp != last) { + if (restart_base_setup_z_and_closest(list, newp)) { + // newp was not dominated by something else. + double new_v = one_contribution_3d(newp); + assert(new_v > 0); + volume += new_v; + add_to_z(newp); + update_links(list, newp); + } + // FIXME: It newp was dominated, can we accumulate the height and update + // hv later? + double height = newp->next[1]->x[3] - newp->x[3]; + assert(height >= 0); + hv += volume * height; + newp = newp->next[1]; + } + return hv; +} + +#endif // _HV4D_PRIV_H diff --git a/c/hv_priv.h b/c/hv_priv.h index 34c7e0b03..320cf870a 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -12,15 +12,25 @@ // ----------------------- Data Structure ------------------------------------- /* - With HV_DIMENSION==3, we have 'struct dlnode * next[1]' instead of 'struct dlnode * next'. - The compiler should be able to remove the extra indirection. */ typedef struct dlnode { - const double * restrict x; // point coordinates (objective vector). + const double * restrict x; // point coordinates (objective vector). +#ifdef HV_RECURSIVE + // In the recursive algorithm, the number of dimensions is not known in + // advance, so next and prev cannot be fixed-size arrays. + struct dlnode ** r_next; // next-node vector for dimensions 5 and above. + struct dlnode ** r_prev; // previous-node vector for dimensions 5 and above. + double * restrict area; // partial area for dimensions 4 and above. + double * restrict vol; // partial volume for dimensions 4 and above. +#endif + /* With HV_DIMENSION==3, we have 'struct dlnode * next[1]' instead of + 'struct dlnode * next'. GCC -O3 is able to remove the extra + indirection, so it is not worth having a special case. */ struct dlnode * next[HV_DIMENSION - 2]; /* keeps the points sorted according to coordinates 2,3 and 4 - (in the case of 2 and 3, only the points swept by 4 are kept) */ + (in the case of 2 and 3, only the points swept by 4 are kept) */ struct dlnode * prev[HV_DIMENSION - 2]; //keeps the points sorted according to coordinates 2 and 3 (except the sentinel 3) + struct dlnode * cnext[2]; //current next #if HV_DIMENSION == 4 || defined(HVC_ONLY) struct dlnode * closest[2]; // closest[0] == cx, closest[1] == cy @@ -31,8 +41,9 @@ typedef struct dlnode { double area, volume; double last_slice_z; // FIXME: Is this really needed? struct dlnode * head[2]; // lowest (x, y) - bool ignore; // hvc should be zero (duplicated or dominated) #endif + // In hvc this is used as a boolean (duplicated or dominated) + dimension_t ignore; // [0, 255] } dlnode_t; // ------------ Update data structure ----------------------------------------- @@ -80,8 +91,9 @@ print_x(const dlnode_t * p) #endif // HV_DIMENSION == 3 // ------------------------ Circular double-linked list ---------------------- + static inline void -reset_sentinels(dlnode_t * restrict list) +reset_sentinels_3d(dlnode_t * restrict list) { dlnode_t * restrict s1 = list; dlnode_t * restrict s2 = list + 1; @@ -89,33 +101,41 @@ reset_sentinels(dlnode_t * restrict list) s1->next[0] = s2; s1->prev[0] = s3; -#if HV_DIMENSION == 4 || defined(HVC_ONLY) - s1->closest[0] = s2; - s1->closest[1] = s1; -#endif -#if HV_DIMENSION == 4 - s1->next[1] = s2; - s1->prev[1] = s3; -#endif s2->next[0] = s3; s2->prev[0] = s1; -#if HV_DIMENSION == 4 || defined(HVC_ONLY) - s2->closest[0] = s2; - s2->closest[1] = s1; -#endif -#if HV_DIMENSION == 4 - s2->next[1] = s3; - s2->prev[1] = s1; -#endif s3->next[0] = s1; s3->prev[0] = s2; + #if HV_DIMENSION == 4 || defined(HVC_ONLY) + s1->closest[0] = s2; + s1->closest[1] = s1; + + s2->closest[0] = s2; + s2->closest[1] = s1; + s3->closest[0] = s2; s3->closest[1] = s1; #endif +} + +static inline void +reset_sentinels(dlnode_t * restrict list) +{ + reset_sentinels_3d(list); + #if HV_DIMENSION == 4 + dlnode_t * restrict s1 = list; + dlnode_t * restrict s2 = list + 1; + dlnode_t * restrict s3 = list + 2; + + s1->next[1] = s2; + s1->prev[1] = s3; + + s2->next[1] = s3; + s2->prev[1] = s1; + s3->next[1] = s1; s3->prev[1] = s2; #endif diff --git a/c/libhv.mk b/c/libhv.mk index fbb0cde17..faa368ab9 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 libmoocore-config.h +LIBHV_HDRS = hv.h hv_priv.h hv3d_priv.h hv4d_priv.h libmoocore-config.h LIBHV_OBJS = $(LIBHV_SRCS:.c=.o) HV_LIB = fpli_hv.a From bdc4689630c888ad7445b74fc67bffd5a8e4bee5 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 14 Nov 2025 11:13:22 +0000 Subject: [PATCH 28/49] DEVEL-README.md, c/README.md: More details about building and testing. --- DEVEL-README.md | 26 ++++++++++++++++++++++++++ c/README.md | 19 ++++++++++++++++--- 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/DEVEL-README.md b/DEVEL-README.md index 4d27893f2..2039211e3 100644 --- a/DEVEL-README.md +++ b/DEVEL-README.md @@ -24,6 +24,32 @@ How to release 1. `git ci -a -m "Start development of v{NEW_VERSION}"` +Building the code +================= + +C code +----- + +See [*Compilation*](https://github.com/multi-objective/moocore/tree/main/c#compilation) to compile the C code under `c/`. + +R and Python packages +--------------------- + +Under `r/` or `python/`, the commands `make build` and `make install` will build and install the R and Python packages. + + +Testing +======== + +The C code and the Python and R packages have their own testsuites. See [C testsuite](https://github.com/multi-objective/moocore/tree/main/c#testsuite) to setup and run the C testsuite. + +The Python and R testsuites are run automatically via Github Actions for every push and pull request. You can run then manually by running `make test` under `r/` or `python/`. + +You can run all testsuites at once by executing at the top-level: + + make test + + How to extend the Python API ============================ diff --git a/c/README.md b/c/README.md index 6e5879ff7..3845bea9f 100644 --- a/c/README.md +++ b/c/README.md @@ -45,9 +45,9 @@ This uses the option `"-march=native"`. If your GCC version does not support `"n make all MARCH=x86-64-v2 -See the [GCC manual](https://gcc.gnu.org/onlinedocs/gcc/Submodel-Options.html** for the names of the architectures supported by your version of GCC. +See the [GCC manual](https://gcc.gnu.org/onlinedocs/gcc/Submodel-Options.html) for the names of the architectures supported by your version of GCC. -You can compile the code with various runtime checks enabled using the option `DEBUG=1`. This will slow down the code significantly. +You can compile the code with various runtime checks enabled using the option `DEBUG=1`. This will slow down the code significantly. If this is too slow for your purposes, you can disable the [sanitizers](https://gcc.gnu.org/onlinedocs/gcc/Instrumentation-Options.html#index-fsanitize_003daddress) using `make all DEBUG=1 SANITIZERS=""`. The build system will try to pick good compiler flags for you, but if you need to, you can override them by passing the option `OPT_CFLAGS`. For example, to disable compiler optimization and enabled debug symbols, you could run: @@ -55,7 +55,7 @@ The build system will try to pick good compiler flags for you, but if you need t If you do not want to see the command line of each compiler -invocation, pass `S=1` to `make`. +invocation, you pass `S=1` to `make`. Embedding (shared library) @@ -314,12 +314,25 @@ The **moocore** executables are validated using a [comprehensive testsuite](http git clone https://github.com/multi-objective/moocore moocore # Download the testsuite git clone https://github.com/multi-objective/testsuite moocore/testsuite +# Install required python packages +python3 -m pip install -r testsuite/requirements.txt # Testing make -C moocore/c/ test # Timing make -C moocore/c/ time ``` +Under `moocore/c/`, `make test` will compile the code with `DEBUG=1 SANTIZERS="" OPT_CFLAGS="-Og -g3"`. This is useful for debugging failures. + +`make time` will compile the code with `DEBUG=0 SANITIZERS="" OPT_CLAGS="-O3 -ftlo -march=native"`. This is useful for benchmarking code and making sure compiler optimizations do not break anything. + +Neither `make test` nor `make time` recompile the code that has not been modified. So you can compile the code first with your preferred compiler flags and then run the testsuite. For example, if you compile with `make all DEBUG=1`, the code will be compiler with several [sanitizers](https://gcc.gnu.org/onlinedocs/gcc/Instrumentation-Options.html#index-fsanitize_003daddress) enabled. You can then call `make test` to run the testsuite with the sanitizers enabled. This is much slower, but very helpful to track crashes. + +Under `c/` there are various targets to run only parts of the testsuite, for example: + + make test-hv + make test-nondominated + ... From eaa51e8bf31355b62e53dfc18d0847867819c44e Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Tue, 25 Nov 2025 17:02:02 +0000 Subject: [PATCH 29/49] Compute a 4D hv contribution in the base case. --- c/hv.c | 35 +++++++++++++++++--- c/hv4d_priv.h | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++ c/hv_priv.h | 3 +- 3 files changed, 123 insertions(+), 6 deletions(-) diff --git a/c/hv.c b/c/hv.c index 0e4abac62..331ce2bf8 100644 --- a/c/hv.c +++ b/c/hv.c @@ -73,6 +73,7 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, head->vol = head->area + d_stop * (n+1); head->x = NULL; // head contains no data head->ignore = 0; // should never get used + head->x_aux = malloc(3*(n+1) * sizeof(double)); // Reserve space for the sentinels. dlnode_t * list4d = new_cdllist(0, ref); @@ -93,6 +94,7 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, 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].x_aux = head->x_aux + i * 3; i++; } } @@ -135,8 +137,12 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } } // Reset x to point to the first objective. - for (i = 0; i < n; i++) + size_t l; + for (i = 0; i < n; i++){ scratch[i]->x -= STOP_DIMENSION; + for(l = 0; l < 3; l++) + scratch[i]->x_aux[l] = scratch[i]->x[l]; + } free(scratch); @@ -152,6 +158,7 @@ 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]); + free(head->x_aux); free_cdllist(head->next[0]); // free 4D sentinels free(head->r_next); free(head->area); @@ -228,6 +235,18 @@ fpli_hv4d(dlnode_t * restrict list, size_t c _attr_maybe_unused) return hv; } +static double +fpli_onec4d(dlnode_t * restrict list, size_t c _attr_maybe_unused, dlnode_t *the_point) +{ + ASSUME(c > 1); + assert(list->next[0] == list->prev[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 contrib = onec4dplusU(list4d, the_point); + return contrib; +} + static double one_point_hv(const double * restrict x, const double * restrict ref, dimension_t d) { @@ -257,7 +276,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); @@ -265,7 +284,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); /* ------------------------------------------------------ @@ -334,7 +353,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. */ @@ -407,7 +429,10 @@ fpli_hv_ge5d(dlnode_t * restrict list, dimension_t dim, size_t c, assert(p1_prev == p1->r_prev[d_stop - 1]); p1->vol[d_stop] = hyperv; assert(p1->ignore == 0); - double hypera = hv_recursive(list, dim - 1, c, ref, bound); + 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. */ diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 0adaddfc0..ec3aaf09a 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -245,4 +245,95 @@ hv4dplusU(dlnode_t * list) return hv; } +static restore_points(dlnode_t * list, dlnode_t * last){ + dlnode_t * newp = (list+1)->next[1]; + while(newp != last){ + for(int i = 0; i < 3; i++){ + newp->x[i] = newp->x_aux[i]; + } + newp = newp->next[1]; + } +} + +/* Compute the hv contribution of "the_point" in d=4 by iteratively computing the one contribution problem in d=3. */ +static double +onec4dplusU(dlnode_t * list, dlnode_t * the_point) +{ + 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]); + + double volume = 0, hv = 0; + dlnode_t * newp = (list+1)->next[1]; + const dlnode_t * last = list+2; + + the_point->closest[0] = list+1; + the_point->closest[1] = list; + + //PART 1: Setup 3D base (TODO: improve) + while (newp != last && newp->x[3] <= the_point->x[3]) { + if(newp != the_point && newp->ignore < 3){ + + if(newp->x[0] <= the_point->x[0] && newp->x[1] <= the_point->x[1] && newp->x[2] <= the_point->x[2]){ + the_point->ignore = 3; + restore_points(list, newp); + return 0; + } + + for(int i = 0; i < 3; i++){ + newp->x[i] = (newp->x[i] >= the_point->x[i]) ? newp->x[i] : the_point->x[i]; + } + + if (restart_base_setup_z_and_closest(list, newp)) { + add_to_z(newp); + update_links(list, newp); + } + } + newp = newp->next[1]; + } + + restart_base_setup_z_and_closest(list, the_point); + double the_point_v = one_contribution_3d(the_point); + assert(the_point_v > 0); + volume = the_point_v; + double height = newp->x[3] - the_point->x[3]; + assert(height >= 0); + hv = volume * height; + + // PART 2: Update the 3D contribution + while (newp != last && + (newp->x[0] > the_point->x[0] || newp->x[1] > the_point->x[1] || newp->x[2] > the_point->x[2])) { + + if (newp != the_point && newp->ignore < 3){ + for(int i = 0; i < 3; i++){ + newp->x[i] = (newp->x[i] >= the_point->x[i]) ? newp->x[i] : the_point->x[i]; + } + if (restart_base_setup_z_and_closest(list, newp)) { + + // newp was not dominated by something else. + double newp_v = one_contribution_3d(newp); + assert(newp_v > 0); + volume -= newp_v; + + add_to_z(newp); + update_links(list, newp); + } + } + // FIXME: It newp was dominated, can we accumulate the height and update + // hv later? + height = newp->next[1]->x[3] - newp->x[3]; + assert(height >= 0); + hv += volume * height; + + newp = newp->next[1]; + } + + restore_points(list, newp); + return hv; +} + #endif // _HV4D_PRIV_H diff --git a/c/hv_priv.h b/c/hv_priv.h index 320cf870a..e54819916 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -15,7 +15,7 @@ */ typedef struct dlnode { - const double * restrict x; // point coordinates (objective vector). + /*const*/ double * restrict x; // point coordinates (objective vector). #ifdef HV_RECURSIVE // In the recursive algorithm, the number of dimensions is not known in // advance, so next and prev cannot be fixed-size arrays. @@ -32,6 +32,7 @@ typedef struct dlnode { struct dlnode * prev[HV_DIMENSION - 2]; //keeps the points sorted according to coordinates 2 and 3 (except the sentinel 3) struct dlnode * cnext[2]; //current next + double * x_aux; // point coordinates (objective vector). #if HV_DIMENSION == 4 || defined(HVC_ONLY) struct dlnode * closest[2]; // closest[0] == cx, closest[1] == cy // FIXME: unused From f11a07fd0bf57fe20312579848df59571326b885 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 26 Nov 2025 09:49:16 +0000 Subject: [PATCH 30/49] * hv4d_priv.h: Add return type. --- c/hv4d_priv.h | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index ec3aaf09a..2a1c0bd62 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -245,10 +245,12 @@ hv4dplusU(dlnode_t * list) return hv; } -static restore_points(dlnode_t * list, dlnode_t * last){ +static inline void +restore_points(dlnode_t * list, dlnode_t * last) +{ dlnode_t * newp = (list+1)->next[1]; - while(newp != last){ - for(int i = 0; i < 3; i++){ + while (newp != last) { + for (int i = 0; i < 3; i++) { newp->x[i] = newp->x_aux[i]; } newp = newp->next[1]; @@ -256,10 +258,10 @@ static restore_points(dlnode_t * list, dlnode_t * last){ } /* Compute the hv contribution of "the_point" in d=4 by iteratively computing the one contribution problem in d=3. */ -static double +static inline double onec4dplusU(dlnode_t * list, dlnode_t * the_point) { - if(the_point->ignore >= 3){ + if (the_point->ignore >= 3) { return 0; } @@ -284,7 +286,7 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) return 0; } - for(int i = 0; i < 3; i++){ + for (int i = 0; i < 3; i++){ newp->x[i] = (newp->x[i] >= the_point->x[i]) ? newp->x[i] : the_point->x[i]; } From 66cfa221a6d676e9d40a52e25bb10894204e3913 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 26 Nov 2025 10:12:07 +0000 Subject: [PATCH 31/49] hv4d_priv.h: Simplify code and add comments. --- c/hv4d_priv.h | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 2a1c0bd62..4918c4357 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -269,25 +269,29 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) assert(list+2 == list->prev[1]); assert(list+1 == list->next[1]); - double volume = 0, hv = 0; dlnode_t * newp = (list+1)->next[1]; const dlnode_t * last = list+2; the_point->closest[0] = list+1; the_point->closest[1] = list; - //PART 1: Setup 3D base (TODO: improve) - while (newp != last && newp->x[3] <= the_point->x[3]) { - if(newp != the_point && newp->ignore < 3){ + const double * the_point_x = the_point->x; + + // PART 1: Setup 3D base (TODO: improve) + while (newp != last && newp->x[3] <= the_point_x[3]) { + // MANUEL: When can newp be equal to the_point? + if (newp != the_point && newp->ignore < 3){ - if(newp->x[0] <= the_point->x[0] && newp->x[1] <= the_point->x[1] && newp->x[2] <= the_point->x[2]){ + if (newp->x[0] <= the_point_x[0] && newp->x[1] <= the_point_x[1] && newp->x[2] <= the_point_x[2]){ the_point->ignore = 3; restore_points(list, newp); return 0; } + // FIXME: This modifies ->x[], why? for (int i = 0; i < 3; i++){ - newp->x[i] = (newp->x[i] >= the_point->x[i]) ? newp->x[i] : the_point->x[i]; + // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; + newp->x[i] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; } if (restart_base_setup_z_and_closest(list, newp)) { @@ -299,20 +303,21 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) } restart_base_setup_z_and_closest(list, the_point); - double the_point_v = one_contribution_3d(the_point); - assert(the_point_v > 0); - volume = the_point_v; - double height = newp->x[3] - the_point->x[3]; + double volume = one_contribution_3d(the_point); + assert(volume > 0); + double height = newp->x[3] - the_point_x[3]; assert(height >= 0); - hv = volume * height; + double hv = volume * height; // PART 2: Update the 3D contribution while (newp != last && - (newp->x[0] > the_point->x[0] || newp->x[1] > the_point->x[1] || newp->x[2] > the_point->x[2])) { + (newp->x[0] > the_point_x[0] || newp->x[1] > the_point_x[1] || newp->x[2] > the_point_x[2])) { - if (newp != the_point && newp->ignore < 3){ + if (newp != the_point && newp->ignore < 3) { + // FIXME: This modifies ->x[], why? for(int i = 0; i < 3; i++){ - newp->x[i] = (newp->x[i] >= the_point->x[i]) ? newp->x[i] : the_point->x[i]; + // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; + newp->x[i] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; } if (restart_base_setup_z_and_closest(list, newp)) { From f00321c84980c621b0d0539e2be94868acc92e22 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Wed, 26 Nov 2025 12:18:28 +0000 Subject: [PATCH 32/49] hv4d_priv.h: Add a few comments. --- c/hv4d_priv.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 4918c4357..87b0467b2 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -261,6 +261,7 @@ restore_points(dlnode_t * list, dlnode_t * last) static inline double onec4dplusU(dlnode_t * list, dlnode_t * the_point) { + // MANUEL: It would be better to move this check to hv.c to avoid calling reset_sentinels_3d. You can add an assert here. if (the_point->ignore >= 3) { return 0; } @@ -313,6 +314,9 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) while (newp != last && (newp->x[0] > the_point_x[0] || newp->x[1] > the_point_x[1] || newp->x[2] > the_point_x[2])) { + // MANUEL: I think newp cannot be equal to the_point here. If it was + // equal, we would have exited the loop. + assert(newp != the_point); // if (newp != the_point && newp->ignore < 3) { // FIXME: This modifies ->x[], why? for(int i = 0; i < 3; i++){ From 96a7aa4b2858320a83232334e8d4e3cb3d85434c Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 28 Nov 2025 18:22:40 +0000 Subject: [PATCH 33/49] hv4d_priv.h: Use weakly_dominates(). --- c/hv4d_priv.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 87b0467b2..01bcab719 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -277,13 +277,14 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) the_point->closest[1] = list; const double * the_point_x = the_point->x; - + // MANUEL: I added this because I think we cannot (and should not) call this function with an empty list. + assert(newp != last); // PART 1: Setup 3D base (TODO: improve) while (newp != last && newp->x[3] <= the_point_x[3]) { // MANUEL: When can newp be equal to the_point? if (newp != the_point && newp->ignore < 3){ - if (newp->x[0] <= the_point_x[0] && newp->x[1] <= the_point_x[1] && newp->x[2] <= the_point_x[2]){ + if (weakly_dominates(newp->x, the_point_x, 3)) { the_point->ignore = 3; restore_points(list, newp); return 0; @@ -307,7 +308,8 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) double volume = one_contribution_3d(the_point); assert(volume > 0); double height = newp->x[3] - the_point_x[3]; - assert(height >= 0); + // It cannot be zero because we exited the loop above. + assert(height > 0); double hv = volume * height; // PART 2: Update the 3D contribution From 7141650caff8498d16717a960b794695fc217f62 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 28 Nov 2025 18:46:33 +0000 Subject: [PATCH 34/49] hv4d_priv.h: Online restore if ignore < 3. --- c/hv4d_priv.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 01bcab719..824b2f059 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -250,8 +250,11 @@ restore_points(dlnode_t * list, dlnode_t * last) { dlnode_t * newp = (list+1)->next[1]; while (newp != last) { - for (int i = 0; i < 3; i++) { - newp->x[i] = newp->x_aux[i]; + // MANUEL: We only modify the points that are not ignored. + if (newp->ignore < 3) { + for (int i = 0; i < 3; i++) { + newp->x[i] = newp->x_aux[i]; + } } newp = newp->next[1]; } From b4dc3a7ec09fe218bbe524732efeaf23807e3764 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Thu, 11 Dec 2025 16:57:14 +0000 Subject: [PATCH 35/49] 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 | 51 +++++++------ c/hv4d_priv.h | 200 ++++++++++++++++++++++++++++++++++++-------------- c/hv_priv.h | 3 +- c/sort.h | 6 ++ 4 files changed, 182 insertions(+), 78 deletions(-) diff --git a/c/hv.c b/c/hv.c index 331ce2bf8..876a43dd0 100644 --- a/c/hv.c +++ b/c/hv.c @@ -73,14 +73,20 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, head->vol = head->area + d_stop * (n+1); head->x = NULL; // head contains no data head->ignore = 0; // should never get used - head->x_aux = malloc(3*(n+1) * sizeof(double)); // Reserve space for the sentinels. dlnode_t * list4d = new_cdllist(0, ref); // 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. + // 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 + dlnode_t * list_aux = (dlnode_t *) malloc((n + 1) * sizeof(*list_aux)); + double * x_aux = malloc(3 * n * sizeof(double)); + list_aux->x = x_aux; + head->prev[0] = list_aux; + list_aux->next[0] = list4d; size_t i = 1; for (size_t j = 0; j < n; j++) { @@ -94,7 +100,6 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, 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].x_aux = head->x_aux + i * 3; i++; } } @@ -106,7 +111,7 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, for (i = 0; i < n; i++) scratch[i] = head + i + 1; - for (int j = d_stop - 2; j >= -1; j--) { + 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 @@ -115,16 +120,17 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, for (i = 0; i < n; i++) scratch[i]->x--; // Sort each dimension independently. + const int j2 = j+2; 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) { + (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]; + scratch[n-1]->next[j2] = list4d+2; + (list4d+2)->prev[j2] = scratch[n-1]; } else { head->r_next[j] = scratch[0]; scratch[0]->r_prev[j] = head; @@ -137,11 +143,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } } // Reset x to point to the first objective. - size_t l; for (i = 0; i < n; i++){ - scratch[i]->x -= STOP_DIMENSION; - for(l = 0; l < 3; l++) - scratch[i]->x_aux[l] = scratch[i]->x[l]; + scratch[i]->x -= STOP_DIMENSION-1; } free(scratch); @@ -157,9 +160,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]); - free(head->x_aux); + assert(head->next[0] == head->prev[0]->next[0]); free_cdllist(head->next[0]); // free 4D sentinels + free(head->prev[0]->x); // free x_aux (4D basecase) + free(head->prev[0]); // free list_aux (4D basecase) free(head->r_next); free(head->area); free(head); @@ -191,6 +195,9 @@ delete_dom(dlnode_t * restrict nodep, dimension_t dim) // Dimension 4. nodep->prev[1]->next[1] = nodep->next[1]; nodep->next[1]->prev[1] = nodep->prev[1]; + // Dimension 3. + nodep->prev[0]->next[0] = nodep->next[0]; + nodep->next[0]->prev[0] = nodep->prev[0]; } static void @@ -214,6 +221,9 @@ reinsert_nobound(dlnode_t * restrict nodep, dimension_t dim) // Dimension 4. nodep->prev[1]->next[1] = nodep; nodep->next[1]->prev[1] = nodep; + // Dimension 3. + nodep->prev[0]->next[0] = nodep; + nodep->next[0]->prev[0] = nodep; } static void @@ -239,11 +249,10 @@ static double fpli_onec4d(dlnode_t * restrict list, size_t c _attr_maybe_unused, dlnode_t *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 contrib = onec4dplusU(list4d, the_point); + dlnode_t * restrict list_aux = list->prev[0]; + double contrib = onec4dplusU(list4d, list_aux, the_point); return contrib; } diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 824b2f059..6405b52cf 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -159,6 +159,93 @@ 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 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], newp->x[3] }; + assert(list+1 == list->next[0]); + dlnode_t * lex_prev = list+1; + dlnode_t * closest1; + 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)); + + 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){ + // check if newp dominates points in the list + assert(p->cnext[0]->x[1] < newx[1]); + assert(p->x[1] >= newx[1]); + while (p->x[0] >= newx[0]) { + remove_from_z(p); + p = p->cnext[1]; + } + // update the closest of points in the list (max. 1), if needed + assert(p->x[0] < newx[0]); + 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) @@ -245,24 +332,11 @@ hv4dplusU(dlnode_t * list) return hv; } -static inline void -restore_points(dlnode_t * list, dlnode_t * last) -{ - dlnode_t * newp = (list+1)->next[1]; - while (newp != last) { - // MANUEL: We only modify the points that are not ignored. - if (newp->ignore < 3) { - for (int i = 0; i < 3; i++) { - newp->x[i] = newp->x_aux[i]; - } - } - newp = newp->next[1]; - } -} + /* 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 * list, dlnode_t * the_point) +onec4dplusU(dlnode_t * list, dlnode_t * list_aux, dlnode_t * the_point) { // MANUEL: It would be better to move this check to hv.c to avoid calling reset_sentinels_3d. You can add an assert here. if (the_point->ignore >= 3) { @@ -272,43 +346,59 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) assert(list+2 == list->prev[0]); assert(list+2 == list->prev[1]); assert(list+1 == list->next[1]); - - dlnode_t * newp = (list+1)->next[1]; + dlnode_t * newp = (list+1)->next[0]; const dlnode_t * last = list+2; - the_point->closest[0] = list+1; - the_point->closest[1] = list; + dlnode_t * z_first = (list+1)->next[0]; + dlnode_t * z_last = last->prev[0]; const double * the_point_x = the_point->x; - // MANUEL: I added this because I think we cannot (and should not) call this function with an empty list. - assert(newp != last); - // PART 1: Setup 3D base (TODO: improve) - while (newp != last && newp->x[3] <= the_point_x[3]) { - // MANUEL: When can newp be equal to the_point? - if (newp != the_point && newp->ignore < 3){ - - if (weakly_dominates(newp->x, the_point_x, 3)) { - the_point->ignore = 3; - restore_points(list, newp); - return 0; - } + double * x_aux = list_aux->x; + size_t c = 0, cx = 0; - // FIXME: This modifies ->x[], why? - for (int i = 0; i < 3; i++){ - // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; - newp->x[i] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; - } + dlnode_t * newp_aux = list_aux+1; //list_aux is a sentinel + + reset_sentinels_3d(list); + restart_list_y(list); - if (restart_base_setup_z_and_closest(list, newp)) { - add_to_z(newp); - update_links(list, newp); + if ((list+1)->next[1] != the_point){ + // PART 1: Setup 3D base (if there are any points below the_point_x[3]) + while (newp != last) { + // MANUEL: When can newp be equal to the_point? + if (newp->x[3] <= the_point_x[3] && newp != the_point && newp->ignore < 3){ + if (newp->x[0] <= the_point_x[0] && newp->x[1] <= the_point_x[1] && newp->x[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; + } + + for (int i = 0; i < 3; i++){ + // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; AG: Yes. x_aux is the coordinate-wise maximum between newp->x and the_point_x + x_aux[cx++] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; + } + + newp_aux->x = x_aux + cx - 3; + if (continue_base_update_z_closest(list, newp_aux, newp_aux->x[2] == the_point_x[2] && c > 0)) { + newp_aux++; + c++; + } } + newp = newp->next[0]; } - newp = newp->next[1]; } + 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]; restart_base_setup_z_and_closest(list, the_point); 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. @@ -319,28 +409,25 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) while (newp != last && (newp->x[0] > the_point_x[0] || newp->x[1] > the_point_x[1] || newp->x[2] > the_point_x[2])) { - // MANUEL: I think newp cannot be equal to the_point here. If it was - // equal, we would have exited the loop. - assert(newp != the_point); // - if (newp != the_point && newp->ignore < 3) { - // FIXME: This modifies ->x[], why? + if (newp->ignore < 3) { for(int i = 0; i < 3; i++){ - // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; - newp->x[i] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; + // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; AG: Yes. x_aux is the coordinate-wise maximum between newp->x and the_point_x + x_aux[cx++] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; } - if (restart_base_setup_z_and_closest(list, newp)) { - + newp_aux->x = x_aux + cx - 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); + double newp_v = one_contribution_3d(newp_aux); assert(newp_v > 0); volume -= newp_v; - add_to_z(newp); - update_links(list, newp); + add_to_z(newp_aux); + update_links(list, newp_aux); } + newp_aux++; } - // FIXME: It newp was dominated, can we accumulate the height and update - // hv later? + // FIXME: If newp was dominated, can we accumulate the height and update + // hv later? AG: Yes height = newp->next[1]->x[3] - newp->x[3]; assert(height >= 0); hv += volume * height; @@ -348,7 +435,10 @@ onec4dplusU(dlnode_t * list, dlnode_t * the_point) newp = newp->next[1]; } - restore_points(list, newp); + // Restore z list + (list+1)->next[0] = z_first; + (list+2)->prev[0] = z_last; + return hv; } diff --git a/c/hv_priv.h b/c/hv_priv.h index e54819916..320cf870a 100644 --- a/c/hv_priv.h +++ b/c/hv_priv.h @@ -15,7 +15,7 @@ */ typedef struct dlnode { - /*const*/ double * restrict x; // point coordinates (objective vector). + const double * restrict x; // point coordinates (objective vector). #ifdef HV_RECURSIVE // In the recursive algorithm, the number of dimensions is not known in // advance, so next and prev cannot be fixed-size arrays. @@ -32,7 +32,6 @@ typedef struct dlnode { struct dlnode * prev[HV_DIMENSION - 2]; //keeps the points sorted according to coordinates 2 and 3 (except the sentinel 3) struct dlnode * cnext[2]; //current next - double * x_aux; // point coordinates (objective vector). #if HV_DIMENSION == 4 || defined(HVC_ONLY) struct dlnode * closest[2]; // closest[0] == cx, closest[1] == cy // FIXME: unused diff --git a/c/sort.h b/c/sort.h index 08488b8f9..f88846988 100644 --- a/c/sort.h +++ b/c/sort.h @@ -44,6 +44,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 025dd6739a0c1e3b5f89fd8b123eb4c849124a5f Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Fri, 12 Dec 2025 11:33:43 +0000 Subject: [PATCH 36/49] *c/hv4d_priv.h: Remove assert and fix minor bug. --- c/hv4d_priv.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 6405b52cf..7f11fba31 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -126,7 +126,9 @@ 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++; - assert(weakly_dominates(px, newx, 4)); + //AG: in onec4dplusU, px may be just a 3-dimensional vector + // so we cannot have this assertion + // assert(weakly_dominates(px, newx, 4)); return false; } @@ -198,16 +200,18 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new // if all points in the list have z-coordinate equal to px[2] if (equalz){ - // check if newp dominates points in the list 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]; } - // update the closest of points in the list (max. 1), if needed + assert(p->x[0] < newx[0]); - p->closest[0] = newp; + // 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); From 3ef5f2d2affe9bc786dc1c687dc40e2669bcf38f Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Fri, 12 Dec 2025 19:21:08 +0000 Subject: [PATCH 37/49] * 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 | 25 +++------- c/hv4d_priv.h | 133 ++++++++++++++++++++++++++++++-------------------- 2 files changed, 86 insertions(+), 72 deletions(-) diff --git a/c/hv.c b/c/hv.c index 876a43dd0..542ebdf6d 100644 --- a/c/hv.c +++ b/c/hv.c @@ -79,12 +79,11 @@ 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 - dlnode_t * list_aux = (dlnode_t *) malloc((n + 1) * sizeof(*list_aux)); - double * x_aux = malloc(3 * n * sizeof(double)); - list_aux->x = x_aux; + // 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; @@ -162,7 +161,7 @@ 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]->x); // free x_aux (4D basecase) + 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); @@ -234,19 +233,7 @@ 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) -{ - ASSUME(c > 1); - assert(list->next[0] == list->prev[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; -} - -static double -fpli_onec4d(dlnode_t * restrict list, size_t c _attr_maybe_unused, dlnode_t *the_point) +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]->next[0]); diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 7f11fba31..d56eb1413 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 @@ -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,9 +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++; - //AG: in onec4dplusU, px may be just a 3-dimensional vector - // so we cannot have this assertion - // assert(weakly_dominates(px, newx, 4)); +#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,20 +169,21 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n } -/* 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). +/** + 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 bool +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], newp->x[3] }; + 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]){ @@ -191,18 +199,18 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new lex_prev = p; // newp->closest[0] = lex_prev - // check if newp is dominated - if (weakly_dominates(p->x, newx, 3)){ + // 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){ + if (equalz) { assert(p->cnext[0]->x[1] < newx[1]); assert(p->x[1] >= newx[1]); - // check if newp dominates points in the list + // Check if newp dominates points in the list. while (p->x[0] >= newx[0]) { remove_from_z(p); p = p->cnext[1]; @@ -220,7 +228,7 @@ continue_base_update_z_closest(dlnode_t * restrict list, dlnode_t * restrict new closest1 = list; - }else{ + } else { //setup z list newp->prev[0] = list->prev[0]->prev[0]; newp->next[0] = list->prev[0]; @@ -307,7 +315,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]); @@ -316,7 +324,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. @@ -337,12 +345,24 @@ hv4dplusU(dlnode_t * list) } +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]); +} -/* Compute the hv contribution of "the_point" in d=4 by iteratively computing the one contribution problem in d=3. */ +// 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 * list, dlnode_t * list_aux, dlnode_t * the_point) +onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, + dlnode_t * restrict the_point) { - // MANUEL: It would be better to move this check to hv.c to avoid calling reset_sentinels_3d. You can add an assert here. + // FIXME: This is never triggered. + assert(the_point->ignore < 3); if (the_point->ignore >= 3) { return 0; } @@ -350,55 +370,58 @@ onec4dplusU(dlnode_t * list, dlnode_t * list_aux, dlnode_t * the_point) 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 * last = list+2; - dlnode_t * z_first = (list+1)->next[0]; + 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]; - const double * the_point_x = the_point->x; - double * x_aux = list_aux->x; - size_t c = 0, cx = 0; - dlnode_t * newp_aux = list_aux+1; //list_aux is a sentinel + 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); - if ((list+1)->next[1] != the_point){ + // 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]) while (newp != last) { - // MANUEL: When can newp be equal to the_point? - if (newp->x[3] <= the_point_x[3] && newp != the_point && newp->ignore < 3){ - if (newp->x[0] <= the_point_x[0] && newp->x[1] <= the_point_x[1] && newp->x[2] <= the_point_x[2]){ + 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) + // Restore modified links (z list). (list+1)->next[0] = z_first; (list+2)->prev[0] = z_last; return 0; } - for (int i = 0; i < 3; i++){ - // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; AG: Yes. x_aux is the coordinate-wise maximum between newp->x and the_point_x - x_aux[cx++] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; - } - - newp_aux->x = x_aux + cx - 3; - if (continue_base_update_z_closest(list, newp_aux, newp_aux->x[2] == the_point_x[2] && c > 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++; - c++; + done_once = true; } } newp = newp->next[0]; } } newp = the_point->next[1]; - while(newp->x[3] <= the_point_x[3]) + 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]; - restart_base_setup_z_and_closest(list, the_point); + // 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; @@ -409,16 +432,17 @@ onec4dplusU(dlnode_t * list, dlnode_t * list_aux, dlnode_t * the_point) assert(height > 0); double hv = volume * height; - // PART 2: Update the 3D contribution - while (newp != last && - (newp->x[0] > the_point_x[0] || newp->x[1] > the_point_x[1] || newp->x[2] > the_point_x[2])) { + // 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) { - for(int i = 0; i < 3; i++){ - // MANUEL: so basically: if (newp->x[i] < the_point_x[i]) newp->x[i] = the_point_x[i]; AG: Yes. x_aux is the coordinate-wise maximum between newp->x and the_point_x - x_aux[cx++] = (newp->x[i] >= the_point_x[i]) ? newp->x[i] : the_point_x[i]; - } - newp_aux->x = x_aux + cx - 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); @@ -428,11 +452,13 @@ onec4dplusU(dlnode_t * list, dlnode_t * list_aux, dlnode_t * the_point) 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] - newp->x[3]; + height = newp->next[1]->x[3] - newpx[3]; assert(height >= 0); hv += volume * height; @@ -445,5 +471,6 @@ onec4dplusU(dlnode_t * list, dlnode_t * list_aux, dlnode_t * the_point) return hv; } +#endif // HV_RECURSIVE #endif // _HV4D_PRIV_H From 5639151aaf0d39656b9739e4e41cf2d4154d88dd Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 13 Dec 2025 12:52:31 +0000 Subject: [PATCH 38/49] Few minor cleanups --- c/hv.c | 63 +++++++++++++++++++++++++++++++-------------------- c/hv4d_priv.h | 5 ++-- 2 files changed, 41 insertions(+), 27 deletions(-) diff --git a/c/hv.c b/c/hv.c index 542ebdf6d..56e2496f5 100644 --- a/c/hv.c +++ b/c/hv.c @@ -181,6 +181,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) { @@ -191,12 +219,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]; - // Dimension 3. - nodep->prev[0]->next[0] = nodep->next[0]; - nodep->next[0]->prev[0] = nodep->prev[0]; + delete_4d(nodep); + delete_3d(nodep); } static void @@ -217,12 +241,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; - // Dimension 3. - nodep->prev[0]->next[0] = nodep; - nodep->next[0]->prev[0] = nodep; + reinsert_4d(nodep); + reinsert_3d(nodep); } static void @@ -413,20 +433,17 @@ 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); + c++; + double hypera = hv_recursive(list, dim - 1, c, ref, bound, p1); - if(dim - 1 == STOP_DIMENSION){ //hypera only has the contribution of 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 @@ -454,10 +471,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 d56eb1413..11a28f8a9 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -389,7 +389,8 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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]) - while (newp != last) { + 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]) { @@ -411,7 +412,7 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, } } newp = newp->next[0]; - } + } while (newp != last); } newp = the_point->next[1]; while (newp->x[3] <= the_point_x[3]) From 15055a7676b062c6c2bf6f6a862b4ac9288d0544 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 39/49] Add documentation. --- c/hv4d_priv.h | 3 ++- python/src/moocore/_moocore.py | 11 ++++++----- r/R/hv.R | 18 +++++++++--------- r/man/hypervolume.Rd | 18 +++++++++--------- 4 files changed, 26 insertions(+), 24 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 11a28f8a9..11f0e8075 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -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) @@ -141,6 +141,7 @@ restart_base_setup_z_and_closest(dlnode_t * restrict list, dlnode_t * restrict n //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]; diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index c25d146de..b5119adb9 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -432,11 +432,12 @@ def hypervolume( for speed. For 5D or higher, it uses a recursive algorithm - :footcite:p:`FonPaqLop06:hypervolume` with HV4D\ :sup:`+` as the base case, - resulting in a :math:`O(n^{d-2})` time complexity and :math:`O(n)` space - complexity in the worst-case, where :math:`d` is the dimension of the - points. Experimental results show that the pruning techniques used may - reduce the time complexity even further. The original proposal + :footcite:p:`FonPaqLop06:hypervolume` that computes 4D contributions + :footcite:p:`GueFon2017hv4d` as the base case, resulting in a + :math:`O(n^{d-2})` time complexity and :math:`O(n)` space complexity in the + worst-case, where :math:`d` is the dimension of the points. 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^{d-2} \log n)`. Andreia P. Guerreiro enhanced the numerical stability of the recursive algorithm by diff --git a/r/R/hv.R b/r/R/hv.R index ae293d35e..bb3f22eb3 100644 --- a/r/R/hv.R +++ b/r/R/hv.R @@ -54,15 +54,15 @@ #' speed. #' #' For 5D or higher, it uses a recursive algorithm -#' \citep{FonPaqLop06:hypervolume} with \eqn{\text{HV4D}^{+}} as the base case, -#' resulting in a \eqn{O(n^{d-2})} time complexity and \eqn{O(n)} space -#' complexity in the worst-case, where \eqn{d} is the dimension of the points. -#' 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^{d-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^{d-2})} +#' time complexity and \eqn{O(n)} space complexity in the worst-case, where +#' \eqn{d} is the dimension of the points. 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^{d-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 870eb6cba..5bffa9fbe 100644 --- a/r/man/hypervolume.Rd +++ b/r/man/hypervolume.Rd @@ -65,15 +65,15 @@ correctly handles weakly dominated points and has been further optimized for speed. For 5D or higher, it uses a recursive algorithm -\citep{FonPaqLop06:hypervolume} with \eqn{\text{HV4D}^{+}} as the base case, -resulting in a \eqn{O(n^{d-2})} time complexity and \eqn{O(n)} space -complexity in the worst-case, where \eqn{d} is the dimension of the points. -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^{d-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^{d-2})} +time complexity and \eqn{O(n)} space complexity in the worst-case, where +\eqn{d} is the dimension of the points. 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^{d-2} \log n)}. Andreia +P. Guerreiro enhanced the numerical stability of the algorithm by avoiding +floating-point comparisons of partial hypervolumes. } \examples{ From 1731299d91277388880eee493a98f83f45816fee 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 40/49] more simplifications --- c/hv.c | 4 ++-- c/hv4d_priv.h | 24 +++++++++--------------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/c/hv.c b/c/hv.c index 56e2496f5..b27364af3 100644 --- a/c/hv.c +++ b/c/hv.c @@ -142,8 +142,8 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, } } // Reset x to point to the first objective. - for (i = 0; i < n; i++){ - scratch[i]->x -= STOP_DIMENSION-1; + for (i = 0; i < n; i++) { + scratch[i]->x -= STOP_DIMENSION - 1; } free(scratch); diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 11f0e8075..2aa9cc922 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -180,11 +180,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]){ @@ -198,13 +196,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] @@ -227,24 +223,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; From 63d2faa9e94b9c82c997f562d9456d624ae843a6 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 13 Dec 2025 16:42:10 +0000 Subject: [PATCH 41/49] avoid newp != the_point. --- c/hv4d_priv.h | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 2aa9cc922..37cc0fb33 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -378,18 +378,20 @@ 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[3] <= the_point_x[3] && 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; + assert(the_point->ignore == 3); // Restore modified links (z list). (list+1)->next[0] = z_first; (list+2)->prev[0] = z_last; @@ -408,6 +410,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]) From 50f8bf87581171a4871a8d26c5ced4454ecbf492 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Sat, 13 Dec 2025 16:51:27 +0000 Subject: [PATCH 42/49] more cleanups --- c/hv4d_priv.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 37cc0fb33..1e2ee6057 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -367,9 +367,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; @@ -390,14 +390,14 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, do { const double * newpx = newp->x; if (newpx[3] <= the_point_x[3] && newp->ignore < 3) { - 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)) { 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; @@ -434,10 +434,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; From 139c27e46329ce395acde684784f9c1a67d26922 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 15 Dec 2025 00:11:53 +0000 Subject: [PATCH 43/49] hv4d_priv.h: Move newp_aux++ within the if(). --- c/hv4d_priv.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 1e2ee6057..1c597b4c3 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -371,7 +371,6 @@ onec4dplusU(dlnode_t * restrict list, dlnode_t * restrict list_aux, 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 @@ -451,10 +450,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 bb2ec51178d62f5a0a429012df4de20feafccb8f Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 15 Dec 2025 00:41:42 +0000 Subject: [PATCH 44/49] hv4d_priv.h: Avoid warning with DEBUG=0. --- c/hv4d_priv.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 1c597b4c3..891a58d04 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -417,9 +417,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; From 3b2089bf9ce22d0a83e1c447952ab931116ac5a9 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Wed, 17 Dec 2025 16:06:47 +0000 Subject: [PATCH 45/49] * 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 891a58d04..1bf381d76 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -233,6 +233,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; @@ -339,6 +340,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) @@ -373,20 +405,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)) { @@ -396,21 +433,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]; @@ -433,7 +513,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 a2f88fa9ceaad06212d26fd2aa7dc954420089b5 Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Thu, 18 Dec 2025 00:10:38 +0000 Subject: [PATCH 46/49] 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 1bf381d76..94b098a8b 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 94559a0388e2712cd28005347254aef342bf94f3 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Mon, 29 Dec 2025 20:04:55 +0000 Subject: [PATCH 47/49] Bump R version to 4.1 --- .github/workflows/R.yml | 3 +-- r/DESCRIPTION | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/R.yml b/.github/workflows/R.yml index 7d2ede166..77e599b17 100644 --- a/.github/workflows/R.yml +++ b/.github/workflows/R.yml @@ -46,8 +46,7 @@ jobs: os: [ ubuntu-22.04, windows-latest, macos-15-intel, macos-14 ] r: [ release ] include: - # Use 4.0 to check with rtools40's older compiler - - { os: windows-latest, r: '4.0' } + - { os: windows-latest, r: '4.1' } # Use latest ubuntu to make it easier to install dependencies - { os: ubuntu-24.04, r: 'devel', http-user-agent: 'release' } env: diff --git a/r/DESCRIPTION b/r/DESCRIPTION index 982a5835b..31bb24d66 100644 --- a/r/DESCRIPTION +++ b/r/DESCRIPTION @@ -16,7 +16,7 @@ Authors@R: c(person("Manuel", "López-Ibáñez", role = c("aut", "cre"), person("Makoto", "Matsumoto", role = "cph", comment = "mt19937 library"), person("Takuji", "Nishimura", role = "cph", comment = "mt19937 library")) Description: Fast implementations of mathematical operations and performance metrics for multi-objective optimization, including filtering and ranking of dominated vectors according to Pareto optimality, hypervolume metric, C.M. Fonseca, L. Paquete, M. López-Ibáñez (2006) , epsilon indicator, inverted generational distance, computation of the empirical attainment function, V.G. da Fonseca, C.M. Fonseca, A.O. Hall (2001) , and Vorob'ev threshold, expectation and deviation, M. Binois, D. Ginsbourger, O. Roustant (2015) , among others. -Depends: R (>= 4.0) +Depends: R (>= 4.1) Imports: matrixStats, Rdpack From 8eff865cb269d7380a053fb462a3d3bc5219b612 Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Tue, 30 Dec 2025 19:46:02 +0000 Subject: [PATCH 48/49] * hv.c: Add FORCE_BUG. --- c/hv.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/c/hv.c b/c/hv.c index b27364af3..815fa780e 100644 --- a/c/hv.c +++ b/c/hv.c @@ -122,6 +122,14 @@ fpli_setup_cdllist(const double * restrict data, dimension_t d, const int j2 = j+2; qsort(scratch, n, sizeof(*scratch), compare_node); if (j <= -1) { +//#define FORCE_BUG +#ifdef FORCE_BUG + if (j2 == 1) { + dlnode_t * tmp = scratch[0]; + scratch[0] = scratch[1]; + scratch[1] = tmp; + } +#endif (list4d+1)->next[j2] = scratch[0]; scratch[0]->prev[j2] = list4d+1; for (i = 1; i < n; i++) { From 11ccef93783519450a5b20a6d613251e5a23f3ee Mon Sep 17 00:00:00 2001 From: "Andreia P. Guerreiro" Date: Wed, 7 Jan 2026 15:09:04 +0000 Subject: [PATCH 49/49] *c/hv4d_priv.h: Fix bug related to repeated coordinates --- c/hv4d_priv.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/c/hv4d_priv.h b/c/hv4d_priv.h index 94b098a8b..1952ad1e2 100644 --- a/c/hv4d_priv.h +++ b/c/hv4d_priv.h @@ -410,10 +410,9 @@ 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] != the_point || the_point->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. @@ -517,7 +516,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.