diff --git a/bibkeys.txt b/bibkeys.txt index 0f901dd9e..7931e1ef6 100644 --- a/bibkeys.txt +++ b/bibkeys.txt @@ -41,6 +41,7 @@ RubMel1998simulation SchEsqLarCoe2012tec SchKer2025r2v2 VelLam1998gp +Vose1991alias WuAza2001metrics ZhoZhaJin2009igdx ZitThi1998ppsn diff --git a/c/NEWS.md b/c/NEWS.md index 389c63de6..0a66bcb2c 100644 --- a/c/NEWS.md +++ b/c/NEWS.md @@ -3,6 +3,10 @@ ## 0.20 * radixsort.h, insort.h: New. + * rng_alias.h: New. Walker-Vose alias method. + * rng.h: New functions for sampling uniformly distributed 32-bits integer in a + bounded interval and for sampling with a given CDF (roulette-wheel method). + * hvapprox.c (hv_approx_fpras): New function. ## 0.19.2 diff --git a/c/common.h b/c/common.h index d1e7deff4..dff5de98c 100644 --- a/c/common.h +++ b/c/common.h @@ -129,4 +129,29 @@ printf_point(const char * prefix, const double * p, dimension_t dim, fprintf(stderr, "%s", suffix); } +static inline double +cumsum_of_vector_double(double * restrict v, size_t n) +{ + ASSUME(n > 0); + for (size_t i = 1; i < n; i++) + v[i] += v[i-1]; + return v[n - 1]; +} + +/** Uses Kahan summation. */ +static inline double +kahan_sum_of_vector_double(const double * restrict v, size_t n) +{ + double sum = 0.0, compensation = 0.0; + for (size_t i = 0; i < n; i++) { + double y = v[i] - compensation; + double z = sum + y; + compensation = (z - sum) - y; + sum = z; + } + return sum; +} + + + #endif /* !MOOCORE_COMMON_H_ */ diff --git a/c/config.h b/c/config.h index e904dd378..3e4046cf0 100644 --- a/c/config.h +++ b/c/config.h @@ -28,6 +28,13 @@ typedef uint_fast8_t dimension_t; #define MOOCORE_HVAPPROX_DIMENSION_MAX 31 #define MOOCORE_HV_DIMENSION_MAX 31 +#ifndef UINT32_MAX +# define UINT32_MAX ((uint32_t)-1) +#endif +#ifndef UINT64_MAX +# define UINT64_MAX ((uint64_t)-1) +#endif + // Use boolvec for boolean arrays that will be passed from/to R/Python. It also // helps with autovectorization. #ifdef R_PACKAGE diff --git a/c/gcc_attribs.h b/c/gcc_attribs.h index 9b22402e0..4f553fa14 100644 --- a/c/gcc_attribs.h +++ b/c/gcc_attribs.h @@ -236,4 +236,15 @@ # define ASAN_POISON_MEMORY_REGION(addr, size) ((void) (addr), (void) (size)) #endif +/* GCC/Clang: use builtins to avoid MinGW isfinite/isnormal macro oddities and + -Werror=float-conversion false positives. +*/ +#if defined(__clang__) || defined(__GNUC__) +# define is_finite(x) __builtin_isfinite(x) +# define is_normal(x) __builtin_isnormal(x) +#else +# define is_finite(x) isfinite(x) // C99 +# define is_normal(x) isnormal(x) // C99 +#endif + #endif /* GCC_ATTRIBUTES */ diff --git a/c/hvapprox.c b/c/hvapprox.c index 6a8636087..b4d7c1666 100644 --- a/c/hvapprox.c +++ b/c/hvapprox.c @@ -5,6 +5,8 @@ #include "hvapprox.h" #include "pow_int.h" #include "rng.h" +#include "rng_alias.h" +#include "sort.h" #define ALMOST_ZERO_WEIGHT 1e-20 @@ -863,3 +865,94 @@ hv_approx_rphi_fang_wang_plus(const double * restrict data, size_t npoints, dime const long double c_m = sphere_area_div_2_pow_d_times_d[dim]; return STATIC_CAST(double, c_m * (expected / STATIC_CAST(long double, nsamples))); } + + +// FIXME: Make sure that this is vectorized. +static double * +compute_vols(const double * restrict points, size_t npoints, dimension_t dim) +{ + double * vol = malloc(npoints * sizeof(*vol)); + for (size_t i = 0; i < npoints; i++) { + const double * restrict pi = points + i * dim; + double tmp = pi[0]; + for (dimension_t d = 1; d < dim; d++) + tmp *= pi[d]; + vol[i] = tmp; + } + return vol; +} + +/** + FPRAS (fully polynomial-time randomized approximation scheme) + + This function returns, with probability :math:`(1 - \delta)`, an + :math:`\epsilon`-approximation of the hypervolume metric with respect to the + given reference point. + + delta : Probability of failure. + + K. Bringmann, T. Friedrich. Approximating the volume of unions and + intersections of high-dimensional geometric objects. Computational Geometry: + Theory and Applications, vol. 43, pages 601-610, 2010. + https://doi.org/10.1016/j.comgeo.2010.03.004 + + This implementation uses Walker-Vose's Alias method to sample from a + discrete distribution in O(1). The naive roulette-wheel approach requires + O(log n). +*/ +double +hv_approx_fpras(const double * restrict data, size_t npoints, dimension_t dim, + const double * restrict ref, const boolvec * restrict maximise, + uint32_t random_seed, double epsilon, double delta) +{ + ASSUME(2 <= dim && dim <= MOOCORE_HVAPPROX_DIMENSION_MAX); + const double * points = transform_and_filter(data, &npoints, dim, ref, maximise); + if (points == NULL) + return 0; + + ASSUME(0 < npoints && npoints < UINT32_MAX); + ASSUME(0 < epsilon && epsilon < 1); + ASSUME(0 < delta && delta < 1); + + const double T_factor = 8 * (log(2) - log(delta)) * (1. + epsilon) / (epsilon * epsilon); + const double T_double = T_factor * (double)npoints; + if (T_double >= (double)UINT64_MAX) + return -1; // This will run for too long! + + const uint64_t T = (uint64_t) ceil(T_double); + double * vols = compute_vols(points, npoints, dim); // VolumeQuery(B_i) + double total_vol = kahan_sum_of_vector_double(vols, npoints); + + rng_state * rng = rng_new(random_seed); + rng_alias_sampler_t * sampler = rng_alias_sampler_new(vols, (uint32_t) npoints); + + uint64_t t_sum = 0, m = 0; + while (true) { + // Sampling should be done in O(1). + uint32_t i = rng_alias_sampler_choose(rng, sampler); + + // SampleQuery(B_i) + double x[MOOCORE_HVAPPROX_DIMENSION_MAX + 1]; + for (dimension_t d = 0; d < dim; d++) + x[d] = rng_random(rng); + const double * restrict p_i = points + i * dim; + for (dimension_t d = 0; d < dim; d++) + x[d] *= p_i[d]; + + while (true) { + if (unlikely(t_sum >= T)) { + free(vols); + rng_free(rng); + rng_alias_sampler_free(sampler); + free((void *) points); + return total_vol * (T_factor / (double) m); + } + t_sum++; + uint32_t j = rng_uniform_u32_ubound(rng, (uint32_t) npoints); + const double * restrict p_j = points + j * dim; + if (weakly_dominates(x, p_j, dim)) // PointQuery(x, B_j) + break; + } + m++; + } +} diff --git a/c/hvapprox.h b/c/hvapprox.h index fd7af7d45..ef41bf2ff 100644 --- a/c/hvapprox.h +++ b/c/hvapprox.h @@ -37,5 +37,18 @@ MOOCORE_API double hv_approx_rphi_fang_wang_plus( const double * restrict ref, const boolvec * restrict maximise, uint_fast32_t nsamples); +/** + FPRAS (fully polynomial-time randomized approximation scheme) + + K. Bringmann, T. Friedrich. Approximating the volume of unions and + intersections of high-dimensional geometric objects. Computational Geometry: + Theory and Applications, vol. 43, pages 601-610, 2010. + https://doi.org/10.1016/j.comgeo.2010.03.004 +*/ +MOOCORE_API double hv_approx_fpras( + const double * restrict data, size_t npoints, dimension_t dim, + const double * restrict ref, const boolvec * restrict maximise, + uint32_t random_seed, double epsilon, double delta); + END_C_DECLS #endif // HV_APPROX_H_ diff --git a/c/main-hvapprox.c b/c/main-hvapprox.c index d59e3a0c1..fde011551 100644 --- a/c/main-hvapprox.c +++ b/c/main-hvapprox.c @@ -4,12 +4,12 @@ --------------------------------------------------------------------- - Copyright (c) 2025 + Copyright (c) 2026 Manuel Lopez-Ibanez 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/. + file, you can obtain one at https://mozilla.org/MPL/2.0/. ---------------------------------------------------------------------- @@ -31,7 +31,7 @@ #include "timer.h" #include "nondominated.h" #include "hvapprox.h" -#define CMDLINE_COPYRIGHT_YEARS "2025" +#define CMDLINE_COPYRIGHT_YEARS "2026" #define CMDLINE_AUTHORS "Manuel Lopez-Ibanez \n" #include "cmdline.h" @@ -39,8 +39,15 @@ static int verbose_flag = 1; static bool union_flag = false; static char *suffix = NULL; -enum approx_method_t { DZ2019_MC=1, DZ2019_HW=2, Rphi_FWEp=3 }; -static const char * approx_method_str[] = {"DZ2019-MC", "DZ2019-HW", "Rphi-FWE+"}; +enum approx_method_t { DZ2019_MC=1, DZ2019_HW=2, Rphi_FWEp=3, FPRAS=4 }; +struct hvapprox_params_t { + enum approx_method_t method; + uint32_t seed; + uint_fast32_t nsamples; + double epsilon; + double delta; +}; +static const char * approx_method_str[] = {"DZ2019-MC", "DZ2019-HW", "Rphi-FWE+", "FPRAS"}; static void usage(void) { @@ -64,12 +71,17 @@ OPTION_VERSION_STR " -s, --suffix=STRING Create an output file for each input file by appending\n" " this suffix. This is ignored when reading from stdin. \n" " If missing, output is sent to stdout. \n" -" -n, --nsamples=N Number of Monte-Carlo samples (N is a positive integer).\n" " -m, --method=M 1: Monte-Carlo sampling using normal distribution; \n" " 2: Hua-Wang deterministic sampling. \n" " 3: Rphi-FWE+ deterministic sampling (default). \n" +" 4: FPRAS (Bringmann, Friedrich, 2010). \n" +" -n, --nsamples=N Number of Monte-Carlo samples (N is a positive integer).\n" OPTION_SEED_STR -" Only method=1. \n" +" Only method=1 or method=4.\n" +" -e, --epsilon=E Desired relative error of the approximation, E > 0.\n" +" Only method=4.\n" +" -d, --delta=D Desired failure probability 0 < D < 1, where (1 - D) \n" +" gives the confidence level. Only method=4.\n" "\n"); } @@ -90,8 +102,7 @@ OPTION_SEED_STR static void hvapprox_file(const char * filename, double * restrict reference, double * restrict maximum, double * restrict minimum, - int * restrict nobj_p, - uint_fast32_t nsamples, enum approx_method_t hv_approx_method, uint32_t seed) + int * restrict nobj_p, struct hvapprox_params_t hvapprox) { double * data = NULL; int * cumsizes = NULL; @@ -150,15 +161,18 @@ hvapprox_file(const char * filename, double * restrict reference, Timer_start (); double volume; - switch (hv_approx_method) { + switch (hvapprox.method) { case DZ2019_MC: - volume = hv_approx_normal(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, nsamples, seed); + volume = hv_approx_normal(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, hvapprox.nsamples, hvapprox.seed); break; case DZ2019_HW: - volume = hv_approx_hua_wang(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, nsamples); + volume = hv_approx_hua_wang(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, hvapprox.nsamples); break; case Rphi_FWEp: - volume = hv_approx_rphi_fang_wang_plus(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, nsamples); + volume = hv_approx_rphi_fang_wang_plus(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, hvapprox.nsamples); + break; + case FPRAS: + volume = hv_approx_fpras(&data[nobj * cumsize], cumsizes[n] - cumsize, nobj, reference, maximise, hvapprox.seed, hvapprox.epsilon, hvapprox.delta); break; default: // LCOV_EXCL_LINE # nocov unreachable(); @@ -186,7 +200,7 @@ hvapprox_file(const char * filename, double * restrict reference, int main(int argc, char *argv[]) { // See the man page for getopt_long for an explanation of these fields. - static const char short_options[] = "hVvqur:s:n:m:S:"; + static const char short_options[] = "hVvqur:s:n:m:S:e:d:"; static const struct option long_options[] = { {"help", no_argument, NULL, 'h'}, {"version", no_argument, NULL, 'V'}, @@ -198,6 +212,8 @@ int main(int argc, char *argv[]) {"method", required_argument, NULL, 'm'}, {"nsamples", required_argument, NULL, 'n'}, {"seed", required_argument, NULL, 'S'}, + {"epsilon", required_argument, NULL, 'e'}, + {"delta", required_argument, NULL, 'd'}, {NULL, 0, NULL, 0} /* marks end of list */ }; @@ -205,9 +221,9 @@ int main(int argc, char *argv[]) double * reference = NULL; int nobj = 0; - uint32_t seed = 0; - uint_fast32_t nsamples = 0; - enum approx_method_t hv_approx_method = Rphi_FWEp; + struct hvapprox_params_t hvapprox = { + .method = Rphi_FWEp, .seed = 0, .nsamples = 0, .epsilon = 0.01, .delta = 0.1 + }; int opt; /* it's actually going to hold a char. */ int longopt_index; @@ -229,33 +245,51 @@ int main(int argc, char *argv[]) case 'n': { // --nsamples char *endp; long int value = strtol(optarg, &endp, 10); - if (endp == optarg || *endp != '\0' || value <= 0 || value == LONG_MAX) { + if (endp == optarg || *endp != '\0' || value <= 0 || value == LONG_MAX) fatal_error("value of --nsamples must be a positive integer '%s'", optarg); - } - nsamples = (uint_fast32_t) value; + hvapprox.nsamples = (uint_fast32_t) value; + break; + } + + case 'e': { // --epsilon + char *endp; + double value = strtod(optarg, &endp); + if (endp == optarg || *endp != '\0' || value <= 0 || value == HUGE_VAL || !is_normal(value)) + fatal_error("value of --epsilon must be a positive floating-point value, not '%s'", optarg); + hvapprox.epsilon = value; + break; + } + + case 'd': { // --delta + char *endp; + double value = strtod(optarg, &endp); + if (endp == optarg || *endp != '\0' || value <= 0 || value >= 1 || !is_normal(value)) + fatal_error("value of --delta must be a floating-point value within (0, 1), not '%s'", optarg); + hvapprox.delta = value; break; } case 'm': // --method switch (*optarg) { case '1': - hv_approx_method = DZ2019_MC; break; + hvapprox.method = DZ2019_MC; break; case '2': - hv_approx_method = DZ2019_HW; break; + hvapprox.method = DZ2019_HW; break; case '3': - hv_approx_method = Rphi_FWEp; break; + hvapprox.method = Rphi_FWEp; break; + case '4': + hvapprox.method = FPRAS; break; default: - fatal_error("valid values of --method (-m) are: 1, 2 or 3, not '%s'", optarg); + fatal_error("valid values of --method (-m) are: 1, 2, 3, or 4 not '%s'", optarg); } break; case 'S': {// --seed char *endp; long int value = strtol(optarg, &endp, 10); - if (endp == optarg || *endp != '\0' || value <= 0) { + if (endp == optarg || *endp != '\0' || value <= 0) fatal_error("value of --seed must be a positive integer '%s'", optarg); - } - seed = (uint32_t) value; + hvapprox.seed = (uint32_t) value; break; } case 'q': // --quiet @@ -271,27 +305,33 @@ int main(int argc, char *argv[]) } } - if (nsamples == 0) + if (hvapprox.nsamples == 0 && hvapprox.method != FPRAS) fatal_error("must specify a value for --nsamples, for example, --nsamples 524288"); + if (hvapprox.nsamples != 0 && hvapprox.method == FPRAS) + fatal_error("--nsamples does not make sense with --method=4"); - if (seed == 0) { - if (hv_approx_method == DZ2019_MC) - seed = (uint32_t) time(NULL); - } else if (hv_approx_method != DZ2019_MC) { - fatal_error("--seed only makes sense with --method=1"); - } + if (hvapprox.method == DZ2019_MC || hvapprox.method == FPRAS) { + if (hvapprox.seed == 0) + hvapprox.seed = (uint32_t) time(NULL); + } else if (hvapprox.seed != 0) + fatal_error("--seed only makes sense with --method=1 or --method=4"); - if (verbose_flag >= 2) - printf("# Method: %s\n# seed: %"PRIu32 "\n# nsamples: %lu\n", - approx_method_str[hv_approx_method - 1], - seed, (unsigned long) nsamples); + if (verbose_flag >= 2) { + printf("# Method: %s\n", approx_method_str[hvapprox.method - 1]); + if (hvapprox.method == DZ2019_MC || hvapprox.method != FPRAS) + printf("# seed: %"PRIu32 "\n", hvapprox.seed); + if (hvapprox.method != FPRAS) + printf("# nsamples: %lu\n", (unsigned long) hvapprox.nsamples); + else + printf("# epsilon: %g\n# delta: %g\n", hvapprox.epsilon, hvapprox.delta); + } int numfiles = argc - optind; if (numfiles < 1) /* Read stdin. */ - hvapprox_file(NULL, reference, NULL, NULL, &nobj, nsamples, hv_approx_method, seed); + hvapprox_file(NULL, reference, NULL, NULL, &nobj, hvapprox); else if (numfiles == 1) { - hvapprox_file (argv[optind], reference, NULL, NULL, &nobj, nsamples, hv_approx_method, seed); + hvapprox_file (argv[optind], reference, NULL, NULL, &nobj, hvapprox); } else { int k; double * maximum = NULL; @@ -312,7 +352,8 @@ int main(int argc, char *argv[]) } } for (k = 0; k < numfiles; k++) - hvapprox_file (argv[optind + k], reference, maximum, minimum, &nobj, nsamples, hv_approx_method, seed); + hvapprox_file (argv[optind + k], reference, maximum, minimum, &nobj, + hvapprox); free(minimum); free(maximum); diff --git a/c/rng.h b/c/rng.h index 00876ec1b..ad1b7bf11 100644 --- a/c/rng.h +++ b/c/rng.h @@ -1,3 +1,6 @@ +#ifndef RNG_H +#define RNG_H + #include "mt19937/mt19937.h" typedef mt19937_state rng_state; @@ -34,9 +37,90 @@ rng_uniform(rng_state * rng, double low, double high) return low + (high - low) * rng_random(rng); } +/** + Returns a uniformly distributed 32-bits integer in [0, n). + + Lemire, Daniel. "Fast Random Integer Generation in an Interval", ACM + Transactions on Modeling and Computer Simulation (TOMACS), 29(1):1-12, + 2019. https://doi.org/10.1145/3230636 + +*/ +static inline uint32_t +rng_uniform_u32_ubound(rng_state * rng, uint32_t n) +{ + uint64_t m = ((uint64_t)mt19937_next32(rng)) * n; + uint32_t leftover = (uint32_t)m; + if (leftover < n) { + // t = 2^32 mod n, expressed without requiring a 64-bit 2^32 value. + const uint32_t t = (uint32_t)(-n) % n; + while (leftover < t) { + m = ((uint64_t)mt19937_next32(rng)) * n; + leftover = (uint32_t)m; + } + } + return (uint32_t)(m >> 32); +} + + +static inline uint32_t +rng_uniform_u32_bounded(rng_state * rng, uint32_t low, uint32_t high) +{ + assert(rng != NULL); + if (low >= high) + return low; + + return low + rng_uniform_u32_ubound(rng, high - low); +} + +static inline void +rng_validate_cdf(const double * cdf _attr_maybe_unused, uint32_t n) +{ + assert(n > 0); + assert(cdf[0] >= 0.0); + for (uint32_t i = 1; i < n; ++i) { + assert(cdf[i - 1] <= cdf[i]); + assert(0.0 <= cdf[i] && cdf[i] <= 1.0); + } + assert(cdf[n - 1] == 1.0); +} + +/** + Returns a random value within [0, n - 1] based on the given CDF using the + roulette-wheel method. + + The CDF values should be non-decreasing and within [0,1]. + + Requires O(log n). +*/ +static inline uint32_t +rng_random_wheel_uint32(rng_state * rng, const double * cdf, uint32_t n) +{ + DEBUG1(rng_validate_cdf(cdf, n)); // Check the CDF is correct. + + double r = rng_random(rng); + if (r < cdf[0]) + return 0; + if (n == 2 || r >= cdf[n - 2]) + return n - 1; + + // Binary search. + uint32_t low = 1, high = n - 2; + while (low < high) { + uint32_t mid = low + (high - low) / 2; + if (r < cdf[mid]) + high = mid; + else + low = mid + 1; + } + + DEBUG1(for (uint32_t j = 0; j < low; j++) assert(cdf[j] < r)); + DEBUG1(for (uint32_t j = low; j < n; j++) assert(r <= cdf[j])); + return low; +} double rng_standard_normal(rng_state *rng); void rng_bivariate_normal_fill(rng_state * rng, double mu1, double mu2, double sigma1, double sigma2, double rho, double *out, int n); +#endif /* RNG_H */ diff --git a/c/rng_alias.h b/c/rng_alias.h new file mode 100644 index 000000000..cbd31c1a3 --- /dev/null +++ b/c/rng_alias.h @@ -0,0 +1,210 @@ +/****************************************************************************** + + O(1)-per-sample discrete distribution sampler built on Walker-Vose's + algorithm for the alias method. + + Given n outcomes with (possibly unnormalized) weights p[0..n-1], the + Walker-Vose algorithm builds a table in O(n) time that then supports drawing + a weighted-random outcome in O(1) time. + + The main interface is: + + rng_alias_sampler_t * + rng_alias_sampler_new(const double * probabilities, uint32_t n); + + void rng_alias_sampler_free(rng_alias_sampler_t * sampler); + + uint32_t rng_alias_sampler_choose(rng_state * rng, + const rng_alias_sampler_t * sampler); + + Based on Keith Schwarz's Java reference implementation: + + https://www.keithschwarz.com/interesting/code/?dir=alias-method + + but with many significant changes. In particular, this implementation uses + two 32-bit random values. It could be even faster if we had a 64-bits RNG. + + See also: + + https://en.wikipedia.org/wiki/Alias_method + https://jugit.fz-juelich.de/mlz/ransampl + https://www.keithschwarz.com/darts-dice-coins/ + +*****************************************************************************/ +#ifndef RNG_ALIAS_SAMPLER_H +#define RNG_ALIAS_SAMPLER_H + +#include +#include +#include "rng.h" + +/** + One row of the alias table. Kept as a single struct (rather than two + parallel arrays) so a sample only ever touches one cache line. +*/ +typedef struct rng_alias_entry { + uint32_t cutoff; + uint32_t alias; // Outcome to return when the coin flip above fails. +} rng_alias_entry_t; + +typedef struct rng_alias_sampler { + uint32_t n; + rng_alias_entry_t table[]; // Flexible array. +} rng_alias_sampler_t; + +#define RNG_ALIAS_RESIDUAL_TOL (64.0 * DBL_EPSILON) + +static inline double +rng_alias_clean_scaled(double q) +{ + if (q < 0.0 && q > -RNG_ALIAS_RESIDUAL_TOL) + return 0.0; + + if (q < 1.0 && q > 1.0 - RNG_ALIAS_RESIDUAL_TOL) + return 1.0; + + if (q > 1.0 && q < 1.0 + RNG_ALIAS_RESIDUAL_TOL) + return 1.0; + + return q; +} + +static inline uint32_t +rng_alias_cutoff(double q) +{ + if (q <= 0.0) + return 0; + + double x = q * 0x1p32; // 2^32 + /* Prevent conversion of an out-of-range floating-point value to uint32_t + if rounding produces a value >= 2^32 - 1. */ + return (x >= (double)UINT32_MAX) ? UINT32_MAX : (uint32_t)x; +} + +static inline bool +size_max_mul_overflows(size_t a, size_t b) +{ + return b != 0 && a > SIZE_MAX / b; +} + +/** + Construct an alias sampler using Vose's algorithm. + + The caller retains ownership of probabilities. + + probabilities do not need to sum up to 1, but must be >= 0. +*/ +static inline rng_alias_sampler_t * +rng_alias_sampler_new(const double * probabilities, uint32_t n) +{ + if (probabilities == NULL || n == 0) + return NULL; + + double sum = kahan_sum_of_vector_double(probabilities, n); + if (sum <= 0 || !is_finite(sum)) + return NULL; + + if (size_max_mul_overflows(n, sizeof(rng_alias_entry_t)) + || size_max_mul_overflows(n, sizeof(double))) { + return NULL; + } + + uint32_t small_end = 0, large_begin = n; + /* Temporary Vose double stack: + [0, small_end) small entries + [large_begin, n) large entries + */ + uint32_t * stack = (uint32_t *) malloc(n * sizeof(*stack)); + double * scaled = (double *) malloc(n * sizeof(*scaled)); + rng_alias_sampler_t * sampler = malloc(sizeof(*sampler) + n * sizeof(rng_alias_entry_t)); + if (stack == NULL || sampler == NULL || scaled == NULL) { + free(stack); + free(scaled); + free(sampler); + return NULL; + } + sampler->n = n; + + for (uint32_t i = 0; i < n; i++) { + double p = probabilities[i]; + if (p < 0 || !is_finite(p)) { + free(stack); + free(scaled); + free(sampler); + return NULL; + } + double q = (p / sum) * (double)n; + q = rng_alias_clean_scaled(q); + scaled[i] = q; + + if (q < 1.0) + stack[small_end++] = i; + else + stack[--large_begin] = i; + } + + while (small_end != 0 && large_begin != n) { + uint32_t small_index = stack[--small_end]; + uint32_t large_index = stack[large_begin++]; + double q = scaled[small_index]; + + sampler->table[small_index].cutoff = rng_alias_cutoff(q); + sampler->table[small_index].alias = large_index; + scaled[large_index] += q - 1.0; + scaled[large_index] = rng_alias_clean_scaled(scaled[large_index]); + + if (scaled[large_index] < 1.0) + stack[small_end++] = large_index; + else + stack[--large_begin] = large_index; + } + + while (large_begin != n) { + uint32_t i = stack[large_begin++]; + sampler->table[i].cutoff = UINT32_MAX; + sampler->table[i].alias = i; + } + + while (small_end != 0) { + uint32_t i = stack[--small_end]; + sampler->table[i].cutoff = UINT32_MAX; + sampler->table[i].alias = i; + } + + free(stack); + free(scaled); + return sampler; +} + + +static inline void +rng_alias_sampler_free(rng_alias_sampler_t * sampler) +{ + free(sampler); +} + + +/** + Sample an integer between [0, sampler->n). + + It consumes two 32-bits random values: + + 1. One for unbiased column selection. + 2. One for the fixed-point alias decision. + + Column selection may very rarely consume another 32-bits value when + rejection is required. + + No floating-point operations are performed. +*/ +static inline uint32_t +rng_alias_sampler_choose(rng_state * rng, const rng_alias_sampler_t * sampler) +{ + assert(rng != NULL && sampler != NULL); + uint32_t column = rng_uniform_u32_ubound(rng, sampler->n); + uint32_t u = mt19937_next32(rng); + const rng_alias_entry_t * entry = sampler->table + column; + return u < entry->cutoff ? column : entry->alias; +} + +#endif /* RNG_ALIAS_SAMPLER_H */ diff --git a/python/benchmarks/bench_hvapprox.py b/python/benchmarks/bench_hvapprox.py index e9537e085..eef24a6e1 100644 --- a/python/benchmarks/bench_hvapprox.py +++ b/python/benchmarks/bench_hvapprox.py @@ -19,25 +19,12 @@ timeit_template_return_1_value, # time_hv_exact ) -from pymoo.indicators.hv.monte_carlo import ( - ApproximateMonteCarloHypervolume as pymoo_hvapprox, -) - # See https://github.com/multi-objective/testsuite/tree/main/data path_to_data = "../../testsuite/data/" assert pathlib.Path(path_to_data).expanduser().exists() files = { - "DTLZLinearShape.3d": dict( - file="DTLZLinearShape.3d.front.1000pts.10", - ref=1, - range=(100, 1000, 100), - ), - "DTLZLinearShape.4d": dict( - file="DTLZLinearShape.4d.front.1000pts.10", - ref=1, - range=(100, 1000, 100), - ), + # It does not make sense to approximate the hypervolume in less than 6D. "DTLZLinearShape.6d": dict( file="DTLZLinearShape.6d.front.700pts.10.xz", ref=1, @@ -106,8 +93,8 @@ def time_hv_exact(name, maxrow, z, ref): "moocore Rphi-FWE+": lambda z, exact: relerror( exact, moocore.hv_approx(z, ref=ref, method="Rphi-FWE+") ), - "pymoo": lambda z, exact, hv=pymoo_hvapprox(ref_point=ref): relerror( - exact, hv.add(z).hv + "moocore FPRAS(eps=0.01, d=0.1)": lambda z, exact: relerror( + exact, moocore.hv_approx_fpras(z, ref=ref) ), } bench = Bench( diff --git a/python/benchmarks/bench_ndom.py b/python/benchmarks/bench_ndom.py index 974b426aa..389e12962 100644 --- a/python/benchmarks/bench_ndom.py +++ b/python/benchmarks/bench_ndom.py @@ -16,9 +16,6 @@ from botorch.utils.multi_objective.pareto import ( is_non_dominated as botorch_is_nondominated, ) -from pymoo.util.nds.non_dominated_sorting import ( - NonDominatedSorting as pymoo_NonDominatedSorting, -) from desdeo.tools.non_dominated_sorting import ( non_dominated as desdeo_is_nondominated, @@ -87,9 +84,6 @@ def get_dataset(name): z, sense=z.shape[1] * ["max"], distinct=False, use_numba=True ) ), - "pymoo": lambda z, nds=pymoo_NonDominatedSorting(): nds.do( - -z, only_non_dominated_front=True - ), "desdeo": lambda z: bool2pos(desdeo_is_nondominated(-z)), "fast_pareto": lambda z: bool2pos( fast_pareto_is_pf(-z, assume_unique_lexsorted=False) @@ -140,11 +134,6 @@ def get_dataset(name): z, sense=z.shape[1] * ["max"], distinct=True, use_numba=True ) ), - # The following packages do not support deduplication so they are - # actually slower because the user needs to remove duplicates. - "pymoo": lambda z, nds=pymoo_NonDominatedSorting(): nds.do( - -z, only_non_dominated_front=True - ), "desdeo": lambda z: bool2pos(desdeo_is_nondominated(-z)), "fast_pareto": lambda z: bool2pos( fast_pareto_is_pf(-z, assume_unique_lexsorted=False) diff --git a/python/doc/source/REFERENCES.bib b/python/doc/source/REFERENCES.bib index de557e2dc..ffe9e72dd 100644 --- a/python/doc/source/REFERENCES.bib +++ b/python/doc/source/REFERENCES.bib @@ -409,6 +409,20 @@ @article{SchKer2025r2v2 doi = {10.1162/evco.a.366} } +@article{Vose1991alias, + author = {Michael D. Vose}, + title = {A linear algorithm for generating random numbers with a given + distribution}, + journal = {IEEE Transactions on Software Engineering}, + year = 1991, + volume = 17, + number = 9, + pages = {972--975}, + annote = {Proposed Walker-Vose alias method for sampling from a + discrete distribution}, + doi = {10.1109/32.92917} +} + @article{WuAza2001metrics, author = {J. Wu and S. Azam}, title = {Metrics for Quality Assessment of a Multiobjective Design diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-time.png deleted file mode 100644 index 69f09df8b..000000000 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-time.png and /dev/null differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-values.png deleted file mode 100644 index 622c211ad..000000000 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.3d-values.png and /dev/null differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-time.png deleted file mode 100644 index 565ae8575..000000000 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-time.png and /dev/null differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-values.png deleted file mode 100644 index ceb97c978..000000000 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.4d-values.png and /dev/null differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png index 71f060d3d..6d014482a 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png index 12b161f85..9dfbb9ef2 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png index d7bc32148..d90c55814 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png index 3460553f0..d42bb31a7 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png index 06de4ea4f..a64307a46 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png index a73683f3f..d9a705cc9 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png index 7f0f853d5..6cec4c3ce 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png index 66e418fe6..133bfdb77 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-ran.6d-time.png b/python/doc/source/_static/bench/hvapprox_bench-ran.6d-time.png index d8c3fd1c2..0673af73b 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-ran.6d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-ran.6d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-ran.6d-values.png b/python/doc/source/_static/bench/hvapprox_bench-ran.6d-values.png index 62c547f33..3831a15fc 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-ran.6d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-ran.6d-values.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-ran.9d-time.png b/python/doc/source/_static/bench/hvapprox_bench-ran.9d-time.png index 2f3281427..172ae8cb0 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-ran.9d-time.png and b/python/doc/source/_static/bench/hvapprox_bench-ran.9d-time.png differ diff --git a/python/doc/source/_static/bench/hvapprox_bench-ran.9d-values.png b/python/doc/source/_static/bench/hvapprox_bench-ran.9d-values.png index 7f51c3c7e..e6d6dbba4 100644 Binary files a/python/doc/source/_static/bench/hvapprox_bench-ran.9d-values.png and b/python/doc/source/_static/bench/hvapprox_bench-ran.9d-values.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-convex-4d-time.png b/python/doc/source/_static/bench/ndom_bench-convex-4d-time.png index a082332b7..9515a606b 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-convex-4d-time.png and b/python/doc/source/_static/bench/ndom_bench-convex-4d-time.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-ran3d-40k-time.png b/python/doc/source/_static/bench/ndom_bench-ran3d-40k-time.png index 61f037850..e87cb1208 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-ran3d-40k-time.png and b/python/doc/source/_static/bench/ndom_bench-ran3d-40k-time.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-rmnk-10d-time.png b/python/doc/source/_static/bench/ndom_bench-rmnk-10d-time.png index ba7deaa87..3ebbfc999 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-rmnk-10d-time.png and b/python/doc/source/_static/bench/ndom_bench-rmnk-10d-time.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-sphere-4d-time.png b/python/doc/source/_static/bench/ndom_bench-sphere-4d-time.png index ff269cfb6..5434a8e87 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-sphere-4d-time.png and b/python/doc/source/_static/bench/ndom_bench-sphere-4d-time.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-sphere-5d-time.png b/python/doc/source/_static/bench/ndom_bench-sphere-5d-time.png index a63487010..12c27558c 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-sphere-5d-time.png and b/python/doc/source/_static/bench/ndom_bench-sphere-5d-time.png differ diff --git a/python/doc/source/_static/bench/ndom_bench-test2D-200k-time.png b/python/doc/source/_static/bench/ndom_bench-test2D-200k-time.png index e31f780f6..38aa98f24 100644 Binary files a/python/doc/source/_static/bench/ndom_bench-test2D-200k-time.png and b/python/doc/source/_static/bench/ndom_bench-test2D-200k-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-10d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-10d-time.png index 3b30813af..a46993906 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-10d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-10d-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-2d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-2d-time.png index 67cb20253..ba6e95bca 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-2d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-2d-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-3d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-3d-time.png index bd788b8fc..cfbcb2ddf 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-3d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-3d-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-4d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-4d-time.png index 02c3aa677..3d08f8a7f 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-4d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-4d-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-5d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-5d-time.png index 797e22299..25ee0f158 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-5d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-5d-time.png differ diff --git a/python/doc/source/_static/bench/ndsort_bench-ran-9d-time.png b/python/doc/source/_static/bench/ndsort_bench-ran-9d-time.png index d408b6994..bbd26a443 100644 Binary files a/python/doc/source/_static/bench/ndsort_bench-ran-9d-time.png and b/python/doc/source/_static/bench/ndsort_bench-ran-9d-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-convex-4d-time.png b/python/doc/source/_static/bench/wndom_bench-convex-4d-time.png index 154c9cdfc..b4992d57f 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-convex-4d-time.png and b/python/doc/source/_static/bench/wndom_bench-convex-4d-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-ran3d-40k-time.png b/python/doc/source/_static/bench/wndom_bench-ran3d-40k-time.png index 9657600ac..c5e8c8255 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-ran3d-40k-time.png and b/python/doc/source/_static/bench/wndom_bench-ran3d-40k-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-rmnk-10d-time.png b/python/doc/source/_static/bench/wndom_bench-rmnk-10d-time.png index 9fd396d44..873b89f4d 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-rmnk-10d-time.png and b/python/doc/source/_static/bench/wndom_bench-rmnk-10d-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-sphere-4d-time.png b/python/doc/source/_static/bench/wndom_bench-sphere-4d-time.png index 3c1405420..df0c1ee2e 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-sphere-4d-time.png and b/python/doc/source/_static/bench/wndom_bench-sphere-4d-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-sphere-5d-time.png b/python/doc/source/_static/bench/wndom_bench-sphere-5d-time.png index f4d3b20b1..03f2cd1ef 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-sphere-5d-time.png and b/python/doc/source/_static/bench/wndom_bench-sphere-5d-time.png differ diff --git a/python/doc/source/_static/bench/wndom_bench-test2D-200k-time.png b/python/doc/source/_static/bench/wndom_bench-test2D-200k-time.png index 757e1f362..e14aaa3d5 100644 Binary files a/python/doc/source/_static/bench/wndom_bench-test2D-200k-time.png and b/python/doc/source/_static/bench/wndom_bench-test2D-200k-time.png differ diff --git a/python/doc/source/benchmarks.inc.rst b/python/doc/source/benchmarks.inc.rst index 2f29e21ad..86c5854d4 100644 --- a/python/doc/source/benchmarks.inc.rst +++ b/python/doc/source/benchmarks.inc.rst @@ -19,10 +19,11 @@ to include in the benchmarks. Not all packages provide the same functionality. For example, `pymoo`_ does not provide the :ref:`epsilon indicator `. `BoTorch`_ and `paretobench`_ only provide the hypervolume. `paretoset`_ and `fast-pareto`_ -only identify nondominated points. `seqme`_ already uses `moocore`_, and -`DESDEO`_, `DEAP`_, `pymoo`_ and `jMetalPy`_ also use `moocore`_ for -hypervolume, but other functionality, such as filtering dominated points, is -still slower than `moocore`_. +only identify nondominated points. `DESDEO`_, `DEAP`_, and `jMetalPy`_ use +`moocore`_ for computing the hypervolume, but other functionality, such as +filtering dominated points, is still slower than `moocore`_. Recent versions +of `seqme`_ and `pymoo`_ (≥0.6.2) already use `moocore`_ for most functionality +benchmarked here. We would like to benchmark `pygmo`_, however, it is currently impossible to install using ``pip`` (See https://github.com/esa/pygmo2/issues/152). @@ -181,17 +182,17 @@ Approximation of the hypervolume -------------------------------- The following plots compare the accuracy and speed of approximating the -hypervolume with the various methods provided by :func:`moocore.hv_approx`. The -plots show that there is no clear winner, in terms of approximation error, -between methods ``Rphi-FWE+`` (default) and ``DZ2019-HW``, but both produce -consistently lower approximation errors than method ``DZ2019-MC`` and than -`pymoo`_. However, ``Rphi-FWE+`` is as fast as ``DZ2019-MC`` and both are consistently faster than ``DZ2019-HW``, in particular with higher number of objectives. The computation time of `pymoo`_ grows rapidly with the number of input points. - -If you compare the plots of **DTLZLinearShape-3d** and **DTLZLinearShape-4d** below to the ones above in the previous section, you can see that the exact computation of the hypervolume in 3D or 4D for thousands of points takes milliseconds, whereas approximating the hypervolume is significantly slower and, thus, not worth doing. - -|hvapprox_bench-DTLZLinearShape-3d-values| |hvapprox_bench-DTLZLinearShape-3d-time| - -|hvapprox_bench-DTLZLinearShape-4d-values| |hvapprox_bench-DTLZLinearShape-4d-time| +hypervolume with the various methods provided by :func:`moocore.hv_approx` and +:func:`moocore.hv_approx_fpras`. The plots show that there is no clear winner, +in terms of approximation error, between methods ``Rphi-FWE+`` (default) and +``DZ2019-HW``, but both produce consistently lower approximation errors than +method ``DZ2019-MC``. However, ``Rphi-FWE+`` is as fast as ``DZ2019-MC`` and +both are consistently faster than ``DZ2019-HW``, in particular with higher +number of objectives. + +The exact computation of the hypervolume for less than 5D takes milliseconds +for thousands of points, whereas approximating the hypervolume is significantly +slower and, thus, not worth doing. Approximating the hypervolume becomes more useful for dimensions higher than 5, where the exact computation becomes noticeably slower with hundreds of points. @@ -210,18 +211,6 @@ For such problems, method ``DZ2019-HW`` becomes significantly slower than |hvapprox_bench-DTLZSphereShape-10d-values| |hvapprox_bench-DTLZSphereShape-10d-time| -.. |hvapprox_bench-DTLZLinearShape-3d-values| image:: _static/bench/hvapprox_bench-DTLZLinearShape.3d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-3d-time| image:: _static/bench/hvapprox_bench-DTLZLinearShape.3d-time.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-4d-values| image:: _static/bench/hvapprox_bench-DTLZLinearShape.4d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-4d-time| image:: _static/bench/hvapprox_bench-DTLZLinearShape.4d-time.png - :width: 49% - .. |hvapprox_bench-DTLZSphereShape-6d-values| image:: _static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png :width: 49% diff --git a/python/doc/source/reference/functions.metrics.rst b/python/doc/source/reference/functions.metrics.rst index b28dc078e..8480ac55a 100644 --- a/python/doc/source/reference/functions.metrics.rst +++ b/python/doc/source/reference/functions.metrics.rst @@ -207,6 +207,7 @@ Approximating the hypervolume metric :toctree: generated/ hv_approx + hv_approx_fpras whv_hype Computing the hypervolume can be time consuming, thus several approaches have diff --git a/python/doc/source/whatsnew/index.rst b/python/doc/source/whatsnew/index.rst index fd3a814ae..6ef09374c 100644 --- a/python/doc/source/whatsnew/index.rst +++ b/python/doc/source/whatsnew/index.rst @@ -10,6 +10,7 @@ Version 0.4.0 - Requires ``numpy>=2.1``. - :func:`~moocore.vorob_t` returns a :class:`~typing.NamedTuple` instead of a dictionary. - :func:`~moocore.is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation. +- :func:`~moocore.hv_approx_fpras` implements the fully polynomial-time randomized approximation scheme by Bringmann and Friedrich. Version 0.3.2 (11/07/2026) diff --git a/python/examples/plot_hv_approx.py b/python/examples/plot_hv_approx.py index 2bf27c447..c15ec6abc 100644 --- a/python/examples/plot_hv_approx.py +++ b/python/examples/plot_hv_approx.py @@ -16,6 +16,7 @@ """ # sphinx_gallery_multi_image = "single" +import time import numpy as np import pandas as pd import matplotlib.pyplot as plt @@ -142,3 +143,62 @@ plt.tight_layout() plt.show() + + +# %% +# +# Fully polynomial-time randomized approximation scheme (FPRAS) +# ------------------------------------------------------------- +# +# :func:`moocore.hv_approx_fpras()` allows obtaining an approximation with +# relative error smaller than :math:`\epsilon` (``epsilon``) with a given +# probability :math:`1-\delta` (``delta``). As the plot below shows, the +# approximation error is often better than the requested value, but the +# computation time increases very quickly for smaller ``epsilon``. +# +shape = "convex-sphere" +ref = 1.1 +npoints = 50 +dim = 6 +reps = 10 +rng = np.random.default_rng(42) +res = [] +for r in range(reps): + z = moocore.generate_ndset(npoints, dim, method=shape, seed=42 + r) + t_start = time.perf_counter() + exact = moocore.hypervolume(z, ref=ref) + t_end = time.perf_counter() - t_start + for epsilon in [0.1, 0.025, 0.01, 0.005]: + for delta in [0.25, 0.1, 0.05]: + t_start = time.perf_counter() + hv = moocore.hv_approx_fpras( + z, ref=ref, seed=rng, epsilon=epsilon, delta=delta + ) + t_end = time.perf_counter() - t_start + res.append( + dict( + r=r, + epsilon=epsilon, + delta=delta, + hverror=np.abs(1 - hv / exact), + time=t_end, + ) + ) + + +df = pd.DataFrame(res) +df["delta"] = df["delta"].astype("category") +df["epsilon"] = df["epsilon"].astype("category") +plt.figure() +ax = sns.boxplot(df, x="epsilon", y="hverror", hue="delta", whis=[0, 100]) +ax.set_title(f"{shape}-{npoints}-{dim}d", fontsize=10) +ax.set(yscale="log", ylabel="Relative error") +ax.yaxis.grid(True) +plt.tight_layout() +plt.figure() +ax = sns.lineplot(data=df, x="epsilon", y="time", hue="delta", marker="o") +ax.set_title(f"{shape}-{npoints}-{dim}d", fontsize=10) +ax.set(yscale="log", ylabel="CPU time (s)") +ax.set_xscale("log", base=10) +plt.tight_layout() +plt.show() diff --git a/python/src/moocore/__init__.py b/python/src/moocore/__init__.py index 393b416b8..b20027e19 100644 --- a/python/src/moocore/__init__.py +++ b/python/src/moocore/__init__.py @@ -15,6 +15,7 @@ filter_dominated_within_sets, generate_ndset, hv_approx, + hv_approx_fpras, hv_contributions, hypervolume, igd, @@ -61,6 +62,7 @@ "get_dataset", "get_dataset_path", "hv_approx", + "hv_approx_fpras", "hv_contributions", "hypervolume", "igd", diff --git a/python/src/moocore/_moocore.py b/python/src/moocore/_moocore.py index f799227ba..00d8e18ce 100644 --- a/python/src/moocore/_moocore.py +++ b/python/src/moocore/_moocore.py @@ -655,7 +655,7 @@ def hypervolume( >>> len(dat) 100 - Dominated points are ignored, so this: + Weakly-dominated points are ignored, so this: >>> moocore.hypervolume(dat, ref=[10, 10]) 93.55331425585321 @@ -998,7 +998,7 @@ def hv_approx( seed: int | np.random.Generator | None = None, method: Literal["DZ2019-HW", "DZ2019-MC", "Rphi-FWE+"] = "Rphi-FWE+", ) -> float: - r"""Approximate the hypervolume indicator. + r"""Approximate the hypervolume indicator via (quasi)-Monte-Carlo sampling. Approximate the value of the hypervolume metric with respect to a given reference point assuming minimization of all objectives. Methods @@ -1100,21 +1100,30 @@ def hv_approx( Merge all the sets of a dataset by removing the set number column: >>> x = moocore.get_dataset("input1.dat")[:, :-1] + >>> len(x) + 100 - Dominated points are ignored, so this: + Weakly-dominated points increase the runtime, but do not change the returned value, + so the following: >>> moocore.hv_approx(x, ref=10) 93.5533 >>> moocore.hv_approx(x, ref=10, method="DZ2019-HW") 93.5533 + >>> moocore.hv_approx(x, ref=10, method="DZ2019-MC", seed=42) + 93.5919 - gives the same hypervolume approximation as this: + give the same hypervolume approximation as: >>> x = moocore.filter_dominated(x) + >>> len(x) + 6 >>> moocore.hv_approx(x, ref=10) 93.5533 >>> moocore.hv_approx(x, ref=10, method="DZ2019-HW") 93.5533 + >>> moocore.hv_approx(x, ref=10, method="DZ2019-MC", seed=42) + 93.5919 The approximation is far from perfect for large number of dimensions: @@ -1146,7 +1155,11 @@ def hv_approx( ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref") - if not is_integer_value(nsamples) or nsamples <= 0 or nsamples > 2147483648: + if ( + not is_integer_value(nsamples) + or nsamples <= 0 + or nsamples >= 2147483648 + ): raise ValueError( f"nsamples ({nsamples}) must be a positive integer value smaller than 2147483648" ) @@ -1172,8 +1185,147 @@ def hv_approx( points_p, npoints, nobj, ref, maximise_p, nsamples ) case _: - raise ValueError("Unknown method = {method}") + raise ValueError(f"Unknown method = {method}") + + return hv + + +@DocSubstitute() +def hv_approx_fpras( + points: ArrayLike, + /, + ref: ArrayLike, + *, + maximise: bool | Sequence[bool] = False, + seed: int | np.random.Generator | None = None, + epsilon: float = 0.01, + delta: float = 0.1, +) -> float: + r"""Approximate the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS). + + This function implements the approximation algorithm by + :footcite:t:`BriFri2010approx`. This algorithm returns, with probability + :math:`(1 - \delta)`, an :math:`\epsilon`-approximation of the hypervolume + metric with respect to the given reference point, assuming minimization of + all objectives by default. + + .. warning:: Lower values of ``epsilon`` (:math:`\epsilon`) or ``delta`` (:math:`\delta`) require significantly longer computation time. + + .. seealso:: For details of the calculation, see the Notes section below. + + See :ref:`Benchmarks: Approximation of the hypervolume `. + + Parameters + ---------- + points : + ${points} + ref : + ${ref_point} + maximise : + ${maximise} + seed : + ${random_seed} + epsilon : + Desired relative error of the approximation, :math:`\epsilon > 0`. + delta : + Desired failure probability :math:`0 < \delta < 1`, :math:`(1 - \delta)` gives the confidence level. + + Returns + ------- + A single numerical value, the approximate hypervolume indicator. + + See Also + -------- + hypervolume, whv_hype, hv_approx + + Notes + ----- + This function computes an approximation :math:`\hat{v}` of the + true hypervolume :math:`v = \text{hyp}_r(A)` of the input points in :math:`A + \subset \mathbb{R}^m` with respect to the reference point :math:`r \in + \mathbb{R}^m`, such that + + .. math:: + \text{Pr}[(1-\epsilon)v \leq \hat{v} \leq (1+\epsilon)v] \geq (1 - \delta) + + where :math:`\epsilon > 0` and :math:`0 < \delta < 1`. + + The algorithm requires :math:`O(\frac{nm}{\epsilon^2}\log\frac{1}{\delta})`. That is, + it is linear on the number of points and dimensions, but quadratic in the + approximation error. In other words, more accurate approximations require + significantly more time. + + In contrast to the (quasi)-Monte-Carlo methods provided by + :func:`~moocore.hv_approx`, the presence of weakly-dominated points not + only increases the runtime, but also changes the returned approximation for + a fixed random seed. + + The implementation uses Walker-Vose's alias method for sampling from a + discrete probability distribution :footcite:p:`Vose1991alias`, which + requires :math:`O(1)` per sample. Using the naive roulette-wheel method + would add, at least, a factor of :math:`O(\log n)` to the above runtime. + + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> x = np.array([[5, 5], [4, 6], [2, 7], [7, 4]]) + >>> moocore.hypervolume(x, ref=[10, 10]) + 38.0 + >>> moocore.hv_approx(x, ref=[10, 10], method="Rphi-FWE+") + 37.99998 + >>> moocore.hv_approx_fpras(x, ref=[10, 10], epsilon=0.1, delta=0.2, seed=42) + 38.1446 + >>> moocore.hv_approx_fpras(x, ref=[10, 10], epsilon=0.01, delta=0.2, seed=42) + 37.9541 + + Contrary to :func:`~moocore.hv_approx`, the presence of dominated points + not only increases the runtime, but also changes the output for the same random + seed, so this: + + >>> x = moocore.get_dataset("input1.dat")[:, :-1] + >>> moocore.hv_approx_fpras(x, ref=10, seed=42) + 93.4812 + + does NOT give the same hypervolume approximation as this: + + >>> x = moocore.filter_dominated(x) + >>> moocore.hv_approx_fpras(x, ref=10, seed=42) + 93.6813 + + """ + # Convert to numpy.array in case the user provides a list. We use + # np.asarray to convert it to floating-point, otherwise if a user inputs + # something like ref = np.array([10, 10]) then numpy would interpret it as + # an int array. + points = np.asarray(points, dtype=float) + nobj = points.shape[1] + ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref") + + if epsilon <= 0: + raise ValueError(f"epsilon must be positive: {epsilon}") + if delta <= 0 or delta >= 1: + raise ValueError(f"delta must be strictly within (0, 1): {delta}") + + maximise_p = _parse_maximise_to_bool_array(maximise, nobj) + points_p, npoints, nobj = np2d_to_double_array( + points, ctype_shape=("size_t", "uint_fast8_t") + ) + ref = ffi.from_buffer("double []", ref) + seed = _get_seed_for_c(seed) + epsilon = ffi.cast("double", float(epsilon)) + delta = ffi.cast("double", float(delta)) + hv = lib.hv_approx_fpras( + points_p, npoints, nobj, ref, maximise_p, seed, epsilon, delta + ) + if hv < 0: + raise ValueError( + f"The requested approximation (epsilon={epsilon}, delta={delta}) would require a very long time" + ) return hv diff --git a/python/src/moocore/libmoocore.h b/python/src/moocore/libmoocore.h index 89b527014..926776d0e 100644 --- a/python/src/moocore/libmoocore.h +++ b/python/src/moocore/libmoocore.h @@ -77,6 +77,10 @@ double hv_approx_rphi_fang_wang_plus(const double * restrict data, const double * restrict ref, const boolvec * restrict maximise, uint_fast32_t nsamples); + +double hv_approx_fpras(const double * restrict data, size_t npoints, dimension_t dim, + const double * restrict ref, const boolvec * restrict maximise, + uint32_t random_seed, double epsilon, double delta); /* typedef ... hype_sample_dist; hype_sample_dist * hype_dist_unif_new(unsigned long seed); diff --git a/python/tests/test_moocore.py b/python/tests/test_moocore.py index 4100ce95b..c69156de8 100644 --- a/python/tests/test_moocore.py +++ b/python/tests/test_moocore.py @@ -553,6 +553,11 @@ def test_hv_approx(dim): signif = 4 if dim < 8 else 3 if dim < 10 else 2 np.testing.assert_approx_equal(true_hv, appr_hv, significant=signif) + appr_hv = moocore.hv_approx_fpras(x, ref=ref) + print(f"{dim}: {(true_hv - appr_hv) / true_hv}") + signif = 4 + np.testing.assert_approx_equal(true_hv, appr_hv, significant=signif) + def test_hv_approx_default_seed(): x = np.full((1, 5), 0.5) @@ -571,13 +576,27 @@ def test_hv_approx_errors(): ValueError, match=r".*must be a positive integer value.*" ): moocore.hv_approx([[0, 0]], [1, 1], method="DZ2019-MC", nsamples="10") - - with pytest.raises(ValueError, match=r".*Unknown method.*"): + with pytest.raises( + ValueError, match=r".*must be a positive integer value smaller than.*" + ): + moocore.hv_approx([[0, 0]], [1, 1], method="DZ2019-MC", nsamples=2**31) + with pytest.raises(ValueError, match=r".*Unknown method.*None"): moocore.hv_approx([[0, 0]], [1, 1], method="None") assert moocore.hv_approx([[1, 1, 1]], ref=1) == 0 +def test_hv_approx_fpras_errors(): + with pytest.raises(ValueError, match=r"epsilon must be positive"): + moocore.hv_approx_fpras([[0, 0]], [1, 1], epsilon=0) + with pytest.raises(ValueError, match=r"delta must be strictly within"): + moocore.hv_approx_fpras([[0, 0]], [1, 1], delta=0) + with pytest.raises(ValueError, match=r"a very long time"): + moocore.hv_approx_fpras([[0, 0]], [1, 1], epsilon=1e-9, delta=0.001) + + assert moocore.hv_approx([[1, 1, 1]], ref=1) == 0 + + def check_hvc(points, ref, err_msg): hvc = moocore.hv_contributions(points, ref=ref) is_nondom = moocore.is_nondominated(points, keep_weakly=True) diff --git a/r/NAMESPACE b/r/NAMESPACE index 8a331f2fa..798f64bd9 100644 --- a/r/NAMESPACE +++ b/r/NAMESPACE @@ -15,6 +15,7 @@ export(epsilon_mult) export(filter_dominated) export(generate_ndset) export(hv_approx) +export(hv_approx_fpras) export(hv_contributions) export(hypervolume) export(igd) diff --git a/r/NEWS.md b/r/NEWS.md index 90caebaad..ace4405b0 100644 --- a/r/NEWS.md +++ b/r/NEWS.md @@ -1,6 +1,7 @@ # moocore (development) * `is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation. + * `hv_approx_fpras()`: New function. Approximates the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS) [@BriFri2010approx]. # moocore 0.3.2 diff --git a/r/R/hv_approx.R b/r/R/hv_approx.R index 820cb4755..7b6eccca2 100644 --- a/r/R/hv_approx.R +++ b/r/R/hv_approx.R @@ -120,3 +120,104 @@ hv_approx <- function(x, reference, maximise = FALSE, nsamples = 262144L, seed = as.integer(nsamples))) } } + + +#' Approximate the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS). +#' +#' This function implements the approximation algorithm by +#' \citet{BriFri2010approx}. This algorithm returns, with probability \eqn{(1 +#' - \delta)}, an \eqn{\epsilon}-approximation of the hypervolume metric with +#' respect to the given reference point, assuming minimization of all +#' objectives by default. Lower values of `epsilon` (\eqn{\epsilon}) or +#' `delta` (\eqn{\delta}) require significantly longer computation time. +#' +#' @details +#' +#' This function computes an approximation \eqn{\hat{v}} of the true +#' hypervolume \eqn{v = \text{hyp}_r(A)} of the input points in \eqn{A +#' \subset \mathbb{R}^m} with respect to the reference point \eqn{r \in +#' \mathbb{R}^m}, such that +#' +#' \deqn{\text{Pr}[(1-\epsilon)v \leq \hat{v} \leq (1+\epsilon)v] \geq (1 - \delta)} +#' +#' where \eqn{\epsilon > 0} and \eqn{0 < \delta < 1}. +#' +#' The algorithm requires \eqn{O(\frac{nm}{\epsilon^2}\log\frac{1}{\delta})}. +#' That is, it is linear on the number of points and dimensions, but quadratic +#' in the approximation error. In other words, more accurate approximations +#' require significantly more time. +#' +#' In contrast to the (quasi)-Monte-Carlo methods provided by [hv_approx()], +#' the presence of weakly-dominated points not only increases the runtime, but +#' also changes the returned approximation for a fixed random seed. +#' +#' The implementation uses Walker-Vose's alias method for sampling from a +#' discrete probability distribution \citep{Vose1991alias}, which requires +#' \eqn{O(1)} per sample. Using the naive roulette-wheel method would add, at +#' least, a factor of \eqn{O(\log n)} to the above runtime. +#' +#' @inheritParams hv_approx +#' +#' @return A single numerical value. +#' +#' @param epsilon `double(1)`\cr Desired relative error of the approximation, +#' \eqn{\epsilon > 0}. +#' +#' @param delta `double(1)`\cr Desired failure probability +#' \eqn{0 < \delta < 1}; \eqn{(1 - \delta)} gives the confidence level. +#' +#' @seealso [hypervolume()], [whv_hype()], [hv_approx()] +#' +#' @author Manuel \enc{López-Ibáñez}{Lopez-Ibanez} +#' +#' @references +#' +#' \insertAllCited{} +#' +#' @doctest +#' +#' x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol=2, byrow=TRUE) +#' @expect equal(38.0) +#' hypervolume(x, ref=10) +#' @expect equal(37.999979) +#' hv_approx(x, ref=10, method="Rphi-FWE+") +#' @expect equal(38.1446, tolerance=1e-4) +#' hv_approx_fpras(x, ref=10, epsilon=0.1, delta=0.2, seed=42) +#' @expect equal(37.9541, tolerance=1e-4) +#' hv_approx_fpras(x, ref=10, epsilon=0.01, delta=0.2, seed=42) +#' +#' @export +#' @concept metrics +hv_approx_fpras <- function(x, reference, maximise = FALSE, seed = NULL, + epsilon = 0.01, delta = 0.1) +{ + x <- as_double_matrix(x) + nobjs <- ncol(x) + + if (!is.numeric(reference)) + stop("a numerical reference vector must be provided") + if (length(reference) == 1L) reference <- rep_len(reference, nobjs) + stopifnot(length(reference) == nobjs) + + if (length(maximise) == 1L) maximise <- rep_len(maximise, nobjs) + stopifnot(length(maximise) == nobjs) + check_dimension_max(nobjs, .libmoocore_constants[["MOOCORE_HVAPPROX_DIMENSION_MAX"]]) + + if (!is.numeric(epsilon) || length(epsilon) != 1L || epsilon <= 0) + stop("epsilon must be a positive numeric value") + if (!is.numeric(delta) || length(delta) != 1L || delta <= 0 || delta >= 1) + stop("delta must be strictly within (0, 1)") + + seed <- if (is.null(seed)) get_seed() else as_integer(seed) + hv <- .Call(hv_approx_fpras_C, + t(x), + as.double(reference), + as.logical(maximise), + seed, + as.double(epsilon), + as.double(delta)) + if (hv < 0) + stop("The requested approximation (epsilon=", epsilon, ", delta=", delta, + ") would require a very long time") + hv +} diff --git a/r/inst/REFERENCES.bib b/r/inst/REFERENCES.bib index de557e2dc..ffe9e72dd 100644 --- a/r/inst/REFERENCES.bib +++ b/r/inst/REFERENCES.bib @@ -409,6 +409,20 @@ @article{SchKer2025r2v2 doi = {10.1162/evco.a.366} } +@article{Vose1991alias, + author = {Michael D. Vose}, + title = {A linear algorithm for generating random numbers with a given + distribution}, + journal = {IEEE Transactions on Software Engineering}, + year = 1991, + volume = 17, + number = 9, + pages = {972--975}, + annote = {Proposed Walker-Vose alias method for sampling from a + discrete distribution}, + doi = {10.1109/32.92917} +} + @article{WuAza2001metrics, author = {J. Wu and S. Azam}, title = {Metrics for Quality Assessment of a Multiobjective Design diff --git a/r/man/hv_approx_fpras.Rd b/r/man/hv_approx_fpras.Rd new file mode 100644 index 000000000..c045d7d62 --- /dev/null +++ b/r/man/hv_approx_fpras.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hv_approx.R +\name{hv_approx_fpras} +\alias{hv_approx_fpras} +\title{Approximate the hypervolume indicator via a fully polynomial-time randomized approximation scheme (FPRAS).} +\usage{ +hv_approx_fpras( + x, + reference, + maximise = FALSE, + seed = NULL, + epsilon = 0.01, + delta = 0.1 +) +} +\arguments{ +\item{x}{\code{matrix()}|\code{data.frame()}\cr Matrix or data frame of numerical +values, where each row gives the coordinates of a point.} + +\item{reference}{\code{numeric()}\cr Reference point as a vector of numerical +values.} + +\item{maximise}{\code{logical()}\cr Whether the objectives must be maximised +instead of minimised. Either a single logical value that applies to all +objectives or a vector of logical values, with one value per objective.} + +\item{seed}{\code{integer(1)}\cr Random seed.} + +\item{epsilon}{\code{double(1)}\cr Desired relative error of the approximation, +\eqn{\epsilon > 0}.} + +\item{delta}{\code{double(1)}\cr Desired failure probability +\eqn{0 < \delta < 1}; \eqn{(1 - \delta)} gives the confidence level.} +} +\value{ +A single numerical value. +} +\description{ +This function implements the approximation algorithm by +\citet{BriFri2010approx}. This algorithm returns, with probability \eqn{(1 +- \delta)}, an \eqn{\epsilon}-approximation of the hypervolume metric with +respect to the given reference point, assuming minimization of all +objectives by default. Lower values of \code{epsilon} (\eqn{\epsilon}) or +\code{delta} (\eqn{\delta}) require significantly longer computation time. +} +\details{ +This function computes an approximation \eqn{\hat{v}} of the true +hypervolume \eqn{v = \text{hyp}_r(A)} of the input points in \eqn{A +\subset \mathbb{R}^m} with respect to the reference point \eqn{r \in +\mathbb{R}^m}, such that + +\deqn{\text{Pr}[(1-\epsilon)v \leq \hat{v} \leq (1+\epsilon)v] \geq (1 - \delta)} + +where \eqn{\epsilon > 0} and \eqn{0 < \delta < 1}. + +The algorithm requires \eqn{O(\frac{nm}{\epsilon^2}\log\frac{1}{\delta})}. +That is, it is linear on the number of points and dimensions, but quadratic +in the approximation error. In other words, more accurate approximations +require significantly more time. + +In contrast to the (quasi)-Monte-Carlo methods provided by \code{\link[=hv_approx]{hv_approx()}}, +the presence of weakly-dominated points not only increases the runtime, but +also changes the returned approximation for a fixed random seed. + +The implementation uses Walker-Vose's alias method for sampling from a +discrete probability distribution \citep{Vose1991alias}, which requires +\eqn{O(1)} per sample. Using the naive roulette-wheel method would add, at +least, a factor of \eqn{O(\log n)} to the above runtime. +} +\references{ +\insertAllCited{} +} +\seealso{ +\code{\link[=hypervolume]{hypervolume()}}, \code{\link[=whv_hype]{whv_hype()}}, \code{\link[=hv_approx]{hv_approx()}} +} +\author{ +Manuel \enc{López-Ibáñez}{Lopez-Ibanez} +} +\concept{metrics} +\examples{ +x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol=2, byrow=TRUE) +hypervolume(x, ref=10) +hv_approx(x, ref=10, method="Rphi-FWE+") +hv_approx_fpras(x, ref=10, epsilon=0.1, delta=0.2, seed=42) +hv_approx_fpras(x, ref=10, epsilon=0.01, delta=0.2, seed=42) +} diff --git a/r/src/Rmoocore.c b/r/src/Rmoocore.c index 78fec90a6..3ed8ae6fb 100644 --- a/r/src/Rmoocore.c +++ b/r/src/Rmoocore.c @@ -512,6 +512,23 @@ hv_approx_rphi_fang_wang_plus_C(SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP N return Rf_ScalarReal(hv); } +SEXP +hv_approx_fpras_C(SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP SEED, SEXP EPSILON, SEXP DELTA) +{ + SEXP_2_DOUBLE_MATRIX(DATA, data, nobj, npoints); + SEXP_2_DOUBLE_VECTOR(REFERENCE, ref, reference_len); + SEXP_2_LOGICAL_INT_VECTOR(MAXIMISE, maximise, maximise_len); + SEXP_2_UINT32(SEED, seed); + + assert(nobj == reference_len); + assert(nobj == maximise_len); + + double epsilon = Rf_asReal(EPSILON); + double delta = Rf_asReal(DELTA); + double hv = hv_approx_fpras(data, npoints, nobj, ref, maximise, seed, epsilon, delta); + return Rf_ScalarReal(hv); +} + #include "r2_exact.h" SEXP diff --git a/r/src/init.h b/r/src/init.h index ec7a6eb8d..f9808d096 100644 --- a/r/src/init.h +++ b/r/src/init.h @@ -21,4 +21,5 @@ DECLARE_CALL(whv_hype_C, SEXP DATA, SEXP IDEAL, SEXP REFERENCE, SEXP NSAMPLES, S DECLARE_CALL(hv_approx_dz2019_mc_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES, SEXP SEED) DECLARE_CALL(hv_approx_dz2019_hw_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES) DECLARE_CALL(hv_approx_rphi_fang_wang_plus_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP NSAMPLES) +DECLARE_CALL(hv_approx_fpras_C, SEXP DATA, SEXP REFERENCE, SEXP MAXIMISE, SEXP SEED, SEXP EPSILON, SEXP DELTA) DECLARE_CALL_VOID(libmoocore_constants) diff --git a/r/tests/testthat/test-doctest-hv_approx_fpras.R b/r/tests/testthat/test-doctest-hv_approx_fpras.R new file mode 100644 index 000000000..cf3a4691f --- /dev/null +++ b/r/tests/testthat/test-doctest-hv_approx_fpras.R @@ -0,0 +1,15 @@ +# Generated by doctest: do not edit by hand +# Please edit file in R/hv_approx.R + +test_that("Doctest: hv_approx_fpras", { + # Created from @doctest for `hv_approx_fpras` + # Source file: R/hv_approx.R + x <- matrix(c(5, 5, 4, 6, 2, 7, 7, 4), ncol = 2, byrow = TRUE) + expect_equal(hypervolume(x, ref = 10), 38) + expect_equal(hv_approx(x, ref = 10, method = "Rphi-FWE+"), 37.999979) + expect_equal(hv_approx_fpras(x, ref = 10, epsilon = 0.1, delta = 0.2, seed = 42), + 38.1446, tolerance = 1e-04) + expect_equal(hv_approx_fpras(x, ref = 10, epsilon = 0.01, delta = 0.2, seed = 42), + 37.9541, tolerance = 1e-04) +}) + diff --git a/r/tests/testthat/test-hv_approx.R b/r/tests/testthat/test-hv_approx.R index ed21e4bc3..197e8115a 100644 --- a/r/tests/testthat/test-hv_approx.R +++ b/r/tests/testthat/test-hv_approx.R @@ -5,6 +5,17 @@ test_that("hv_approx errors", { expect_equal(hv_approx(x, ref = 1), 0) }) +test_that("hv_approx_fpras errors", { + x <- matrix(c(0, 0), ncol = 2) + expect_error(hv_approx_fpras(x, reference = c(1, 1), epsilon = 0), + "epsilon must be a positive numeric value") + expect_error(hv_approx_fpras(x, reference = c(1, 1), delta = 0), + "delta must be strictly within") + expect_error(hv_approx_fpras(x, reference = c(1, 1), epsilon = 1e-9, delta = 0.001), + "would require a very long time") + expect_equal(hv_approx_fpras(matrix(c(2, 2), ncol = 2), reference = c(1, 1)), 0) +}) + for (dim in seq(3L, 10L)) { test_that(paste0("hv_approx dim=", dim), { x <- matrix(replicate(dim, 0.5), ncol=dim) @@ -27,5 +38,11 @@ for (dim in seq(3L, 10L)) { expect_equal(hv_approx(x, ref=ref), true_hv, tolerance = 10**-signif, info = paste0("dim=", dim, ", signif=", signif, " error=", log10(abs(true_hv - appr_hv) / true_hv))) + + appr_hv <- hv_approx_fpras(x, ref=ref) + signif <- 4 + expect_equal(appr_hv, true_hv, + tolerance = 10**-signif, info = paste0("dim=", dim, ", signif=", signif, + " error=", log10(abs(true_hv - appr_hv) / true_hv))) }) } diff --git a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.4d.rds b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.4d.rds index 611d3a3ce..92baac2dd 100644 Binary files a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.4d.rds and b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.4d.rds differ diff --git a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.5d.rds b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.5d.rds index a64dd10f2..2fea6a763 100644 Binary files a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.5d.rds and b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.5d.rds differ diff --git a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.6d.rds b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.6d.rds index f2e5d2bd3..3ec0118b1 100644 Binary files a/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.6d.rds and b/r/vignettes/articles/bench/bench-hv-DTLZLinearShape.6d.rds differ diff --git a/r/vignettes/articles/bench/bench-hv-spherical-3d.rds b/r/vignettes/articles/bench/bench-hv-spherical-3d.rds index bdc66eb9a..064e7d154 100644 Binary files a/r/vignettes/articles/bench/bench-hv-spherical-3d.rds and b/r/vignettes/articles/bench/bench-hv-spherical-3d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-convex-3d.rds b/r/vignettes/articles/bench/bench-ndom-convex-3d.rds index 6ab40d3a9..824f4e69e 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-convex-3d.rds and b/r/vignettes/articles/bench/bench-ndom-convex-3d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-convex-4d.rds b/r/vignettes/articles/bench/bench-ndom-convex-4d.rds index 40c9c24c3..3cd06a4e8 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-convex-4d.rds and b/r/vignettes/articles/bench/bench-ndom-convex-4d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-ran3d-40k.rds b/r/vignettes/articles/bench/bench-ndom-ran3d-40k.rds index 753d33dc2..3cdd2766b 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-ran3d-40k.rds and b/r/vignettes/articles/bench/bench-ndom-ran3d-40k.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-ran4d.rds b/r/vignettes/articles/bench/bench-ndom-ran4d.rds index b4b682610..abaee9f87 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-ran4d.rds and b/r/vignettes/articles/bench/bench-ndom-ran4d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-rmnk-10d.rds b/r/vignettes/articles/bench/bench-ndom-rmnk-10d.rds index 826811e08..65ffa8bef 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-rmnk-10d.rds and b/r/vignettes/articles/bench/bench-ndom-rmnk-10d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-sphere-3d.rds b/r/vignettes/articles/bench/bench-ndom-sphere-3d.rds index 76e1d3d66..37fde7271 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-sphere-3d.rds and b/r/vignettes/articles/bench/bench-ndom-sphere-3d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-sphere-5d.rds b/r/vignettes/articles/bench/bench-ndom-sphere-5d.rds index fafdad47d..e599ac574 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-sphere-5d.rds and b/r/vignettes/articles/bench/bench-ndom-sphere-5d.rds differ diff --git a/r/vignettes/articles/bench/bench-ndom-test2D-200k.rds b/r/vignettes/articles/bench/bench-ndom-test2D-200k.rds index 622a60d26..726c1ecf9 100644 Binary files a/r/vignettes/articles/bench/bench-ndom-test2D-200k.rds and b/r/vignettes/articles/bench/bench-ndom-test2D-200k.rds differ diff --git a/r/vignettes/articles/benchmarks.Rmd b/r/vignettes/articles/benchmarks.Rmd index 12195257b..9fc6ffd8f 100644 --- a/r/vignettes/articles/benchmarks.Rmd +++ b/r/vignettes/articles/benchmarks.Rmd @@ -30,7 +30,7 @@ knitr::opts_chunk$set( run_benchmarks <- params$run_benchmarks ``` -The following plots compare the performance of [`moocore`][moocore] against [`emoa`][emoa] and [`bbotk`][bbotk]. Other R packages are not included in the comparison because they are based on these packages for the functionality benchmarked, so they are **at least as slow** as them. For example [`GPareto`][GPareto], [`mlr3mbo`][mlr3mbo], [`rmoo`][rmoo] and [`bbotk`][bbotk] use [`emoa`][emoa] to compute the hypervolume. Not all packages provide the same functionality. +The following plots compare the performance of [`moocore`][moocore] against [`emoa`][emoa], [`targeted`][targeted] and [`caRamel`][caRamel]. Other R packages are not included in the comparison because they are based on these packages for the functionality benchmarked, so they are **at least as slow** as them. For example, [`GPareto`][GPareto] and [`rmoo`][rmoo] use [`emoa`][emoa] to compute the hypervolume. [`bbotk`][bbotk] already uses [`moocore`][moocore]. Not all packages provide the same functionality.
Show benchmarking setup code @@ -117,17 +117,23 @@ benchmark_plot <- function (x, title = "", only_seconds=TRUE, ...) { get_package_version <- function(package) paste0(package, " (", as.character(packageVersion(package)), ")") -benchmark <- function(name, x, N, setup, expr.list, prefix, title) { +benchmark <- function(name, x, N, setup, expr.list, prefix, title, check) +{ rds_file <- paste0("bench/bench-", prefix, "-", name, ".rds") if (run_benchmarks || !file.exists(rds_file)) { lapply(names(expr.list), library, character.only = TRUE) names(expr.list) <- sapply(names(expr.list), get_package_version, USE.NAMES=FALSE) + if (!requireNamespace("atime", quietly = TRUE) + || packageVersion("atime") < "2026.7.1") { + stop("Please install atime >= 2026.7.1") + } res <- substitute(atime::atime( N = N, expr.list = expr.list, setup = SETUP, result=FALSE, times=5, + check=check, seconds.limit=10), list(SETUP=setup)) res <- eval(res) saveRDS(res, file = rds_file) @@ -144,12 +150,12 @@ get_ndset <- function(x, filter) do.call(moocore::generate_ndset, x$generate) } -benchmark_all <- function(files, prefix, title, setup, expr.list, filter) +benchmark_all <- function(files, prefix, title, setup, expr.list, filter, check=TRUE) { for (name in names(files)) { p <- benchmark(name = name, x = get_ndset(files[[name]], filter=filter), N = files[[name]]$N, prefix = prefix, title = title, - setup = setup, expr.list = expr.list) + setup = setup, expr.list = expr.list, check=check) print(p) } } @@ -161,7 +167,7 @@ benchmark_all <- function(files, prefix, title, setup, expr.list, filter) The following plots compare the speed of finding (non)dominated solutions, equivalent to `moocore::is_nondominated()`, in 2D, 3D, 4D, 5D and 10D. The plots show that [`moocore`][moocore] is always 10 times faster than [`GPGame`][GPGame], -[`targeted`][targeted] and [`bbotk`][bbotk] for any number of objectives. In addition, +[`caRamel`][caRamel] and [`targeted`][targeted] for any number of objectives. In addition, [`GPGame`][GPGame] calculates wrong values with repeated coordinates ([vpicheny/GPGame#2](https://github.com/vpicheny/GPGame/issues/2)). The `rand4d` testcase is somewhat different from the others, as it is mostly composed of dominated points (only 10590 points are @@ -184,7 +190,6 @@ setup <- quote({ expr.list <- list( moocore = quote(which(moocore::is_nondominated(z))), - bbotk = quote(which(bbotk::is_dominated(tz))), GPGame = quote(GPGame::nonDom(z)), targeted = quote(targeted::nondom(z)), caRamel = quote(caRamel::pareto(nz)) @@ -202,7 +207,7 @@ files <- list( ) benchmark_all(files, prefix="ndom", title = "is_(non)dominated()", - setup = setup, expr.list = expr.list, filter=FALSE) + setup = setup, expr.list = expr.list, filter=FALSE, check=FALSE) ``` @@ -248,7 +253,7 @@ benchmark_all(files, prefix="hv", title = "HV Computation", setup = setup, expr.list = expr.list, filter=TRUE) ``` -As the plots show, [`moocore`][moocore] is always faster than [`emoa`][emoa] and, hence, faster than [`GPareto`][GPareto], [`mlr3mbo`][mlr3mbo], [`rmoo`][rmoo] and [`bbotk`][bbotk]. +As the plots show, [`moocore`][moocore] is always faster than [`emoa`][emoa] and, hence, faster than [`GPareto`][GPareto] and [`rmoo`][rmoo]. # Hypervolume contribution @@ -258,6 +263,7 @@ The only R package, other than [`moocore`][moocore], able to compute hypervolume [GPGame]: https://cran.r-project.org/package=GPGame [GPareto]: http://cran.r-project.org/package=GPareto [bbotk]: https://cran.r-project.org/package=bbotk +[caRamel]: https://cran.r-project.org/package=caRamel [emoa]: https://cran.r-project.org/package=emoa [mlr3mbo]: https://cran.r-project.org/package=mlr3mbo [rmoo]: https://cran.r-project.org/package=rmoo diff --git a/r/vignettes/articles/hv_approx.Rmd b/r/vignettes/articles/hv_approx.Rmd index a46c8b73c..458f7da59 100644 --- a/r/vignettes/articles/hv_approx.Rmd +++ b/r/vignettes/articles/hv_approx.Rmd @@ -158,4 +158,60 @@ ggplot(df, aes(x = samples, y = hverror, color = Method)) + theme_bw() ``` +# FPRAS: Fully Polynomial-Time Randomized Approximation Scheme + +`hv_approx_fpras()` allows obtaining an approximation with relative error +smaller than $\epsilon$ (`epsilon`) with a given probability $1-\delta$ +(`delta`). As the plot below shows, the approximation error is often better +than the requested value, but the computation time increases very quickly for +smaller `epsilon`. + +```{r fpras} +#| fig.alt = "Boxplot of relative error and line plot of CPU time for FPRAS with varying epsilon and delta." +library(ggplot2) +library(scales) +shape <- "convex-sphere" +ref <- 1.1 +npoints <- 50 +dim <- 6 +nreps <- 10 +set.seed(42) +res <- NULL +for (r in seq_len(nreps)) { + z <- generate_ndset(npoints, dim, method = shape, seed = 42 + r) + exact <- hypervolume(z, reference = ref) + for (epsilon in c(0.1, 0.025, 0.01, 0.005)) { + for (delta in c(0.25, 0.1, 0.05)) { + t_start <- proc.time()[["elapsed"]] + hv <- hv_approx_fpras(z, reference = ref, epsilon = epsilon, delta = delta) + t_end <- proc.time()[["elapsed"]] - t_start + res <- rbind(res, data.frame( + r = r, + epsilon = epsilon, + delta = delta, + hverror = abs(1 - hv / exact), + time = t_end + )) + } + } +} +res[["delta"]] <- as.factor(res[["delta"]]) + +ggplot(res, aes(x = as.factor(epsilon), y = hverror, fill = delta)) + + geom_boxplot() + + scale_y_log10(labels = label_log()) + + labs(x = "epsilon", y = "Relative error") + + ggtitle(paste0(shape, "-", npoints, "-", dim, "d")) + + theme_bw() + +ggplot(res, aes(x = epsilon, y = time, color = delta, group = delta)) + + stat_summary(fun = mean, geom = "line") + + stat_summary(fun = mean, geom = "point") + + scale_y_log10(labels = label_log()) + + scale_x_log10(breaks = sort(unique(res[["epsilon"]]))) + + labs(y = "CPU time (s)") + + ggtitle(paste0(shape, "-", npoints, "-", dim, "d")) + + theme_bw() +``` + # References diff --git a/tools/test-rng_alias.c b/tools/test-rng_alias.c new file mode 100644 index 000000000..de0301c85 --- /dev/null +++ b/tools/test-rng_alias.c @@ -0,0 +1,909 @@ + +/* + * test-harness.c + * Statistical and correctness tests for a Walker-Vose alias sampler. + * The random-wheel reference is assumed to be provided by rng.h. + */ +#include +#include +#include +#include +#include +#include + +#include "rng.h" +#include "rng_alias.h" + +#define MASTER_SEED UINT32_C(0x5EED1234) +#define STANDARD_SAMPLES UINT64_C(1000000) +#define RANDOM_SAMPLES UINT64_C(250000) +#define SERIAL_SAMPLES UINT64_C(500000) +#define LARGE_TABLE_SAMPLES UINT64_C(750000) +#define LIFECYCLE_ITERATIONS UINT32_C(20000) +#define RANDOM_CASES UINT32_C(500) +#define RANDOM_MAX_N UINT32_C(257) +#define STAT_P_MIN 1.0e-4 +#define RANDOM_INDIVIDUAL_MIN 1.0e-10 +#define FISHER_P_MIN 1.0e-4 +#define MIN_EXPECTED 5.0 + +static unsigned g_tests_run; +static unsigned g_tests_failed; + +struct fisher_accumulator { + double statistic; + unsigned count; +}; + +struct sample_result { + uint64_t *counts; + uint64_t samples; + uint64_t out_of_range; + uint64_t zero_hits; +}; + +struct gof_result { + double pearson; + double pearson_p; + double g; + double g_p; + int df; +}; + +struct two_sample_result { + double chi2; + double pvalue; + double max_diff; + int df; +}; + +static void record_result(const char *name, int ok) +{ + ++g_tests_run; + if (!ok) + ++g_tests_failed; + printf("%-42s %s\n", name, ok ? "PASS" : "FAIL"); +} + +static int is_weight_vector(const double *w, uint32_t n) +{ + double sum = 0.0; + uint32_t i; + + if (w == NULL || n == 0) + return 0; + for (i = 0; i < n; ++i) { + if (!isfinite(w[i]) || w[i] < 0.0) + return 0; + sum += w[i]; + } + return isfinite(sum) && sum > 0.0; +} + +static int normalize_weights(const double *w, double *p, uint32_t n) +{ + double sum = 0.0; + uint32_t i; + + if (!is_weight_vector(w, n) || p == NULL) + return 0; + for (i = 0; i < n; ++i) + sum += w[i]; + if (!(sum > 0.0) || !isfinite(sum)) + return 0; + for (i = 0; i < n; ++i) + p[i] = w[i] / sum; + return 1; +} + +static int make_cdf(const double *w, double *cdf, uint32_t n) +{ + double *p; + double c = 0.0; + uint32_t i; + + if (!is_weight_vector(w, n) || cdf == NULL) + return 0; + p = (double *)malloc((size_t)n * sizeof(*p)); + if (p == NULL) + return 0; + if (!normalize_weights(w, p, n)) { + free(p); + return 0; + } + for (i = 0; i < n; ++i) { + c += p[i]; + cdf[i] = c; + } + cdf[n - 1] = 1.0; + free(p); + return 1; +} + +static uint64_t sum_counts(const uint64_t *counts, uint32_t n) +{ + uint64_t total = 0; + uint32_t i; + + for (i = 0; i < n; ++i) + total += counts[i]; + return total; +} + +static double u32_to_open01(uint32_t x) +{ + return ((double)x + 0.5) / 4294967296.0; +} + +static double gamma_q(double a, double x) +{ + const int max_iter = 10000; + const double eps = 3.0e-14; + const double fpmin = DBL_MIN / eps; + int i; + + if (a <= 0.0 || x < 0.0) + return NAN; + if (x == 0.0) + return 1.0; + if (x < a + 1.0) { + double ap = a; + double sum = 1.0 / a; + double del = sum; + + for (i = 1; i <= max_iter; ++i) { + ++ap; + del *= x / ap; + sum += del; + if (fabs(del) < fabs(sum) * eps) + break; + } + return 1.0 - sum * exp(-x + a * log(x) - lgamma(a)); + } + else { + double b = x + 1.0 - a; + double c = 1.0 / fpmin; + double d = 1.0 / b; + double h = d; + + for (i = 1; i <= max_iter; ++i) { + double an = -(double)i * ((double)i - a); + double del; + + b += 2.0; + d = an * d + b; + if (fabs(d) < fpmin) + d = fpmin; + c = b + an / c; + if (fabs(c) < fpmin) + c = fpmin; + d = 1.0 / d; + del = d * c; + h *= del; + if (fabs(del - 1.0) < eps) + break; + } + return exp(-x + a * log(x) - lgamma(a)) * h; + } +} + +static double chi_square_pvalue(double stat, int df) +{ + if (df <= 0 || !isfinite(stat)) + return NAN; + + /* + * Pearson and G statistics are non-negative mathematically, but the + * accumulated G statistic can end up as a tiny negative number from + * floating-point roundoff when observed and expected counts are very + * close. Treat that as exact zero rather than as an invalid statistic. + */ + if (stat < 0.0) { + if (stat > -1.0e-10) + stat = 0.0; + else + return NAN; + } + + return gamma_q(0.5 * (double)df, 0.5 * stat); +} + +static void fisher_init(struct fisher_accumulator *f) +{ + f->statistic = 0.0; + f->count = 0; +} + +static void fisher_add(struct fisher_accumulator *f, double p) +{ + if (!isfinite(p) || p <= 0.0) + p = DBL_MIN; + if (p > 1.0) + p = 1.0; + f->statistic += -2.0 * log(p); + ++f->count; +} + +static double fisher_pvalue(const struct fisher_accumulator *f) +{ + if (f->count == 0) + return NAN; + return chi_square_pvalue(f->statistic, 2 * (int)f->count); +} + +static int alloc_sample_result(struct sample_result *r, uint32_t n) +{ + r->counts = (uint64_t *)calloc(n, sizeof(*r->counts)); + r->samples = 0; + r->out_of_range = 0; + r->zero_hits = 0; + return r->counts != NULL; +} + +static void free_sample_result(struct sample_result *r) +{ + free(r->counts); + r->counts = NULL; +} + +static int sample_vose(rng_state *rng, rng_alias_sampler_t *sampler, + const double *p, uint32_t n, + uint64_t samples, struct sample_result *r) +{ + uint64_t i; + + memset(r->counts, 0, (size_t)n * sizeof(*r->counts)); + r->samples = samples; + r->out_of_range = 0; + r->zero_hits = 0; + for (i = 0; i < samples; ++i) { + uint32_t x = rng_alias_sampler_choose(rng, sampler); + + if (x >= n) { + ++r->out_of_range; + continue; + } + if (p[x] == 0.0) + ++r->zero_hits; + ++r->counts[x]; + } + return r->out_of_range == 0 && sum_counts(r->counts, n) == samples; +} + +static int sample_wheel(rng_state *rng, const double *w, const double *p, + uint32_t n, uint64_t samples, + struct sample_result *r) +{ + double *cdf; + uint64_t i; + + cdf = (double *)malloc((size_t)n * sizeof(*cdf)); + if (cdf == NULL) + return 0; + if (!make_cdf(w, cdf, n)) { + free(cdf); + return 0; + } + memset(r->counts, 0, (size_t)n * sizeof(*r->counts)); + r->samples = samples; + r->out_of_range = 0; + r->zero_hits = 0; + for (i = 0; i < samples; ++i) { + uint32_t x = rng_random_wheel_uint32(rng, cdf, n); + + if (x >= n) { + ++r->out_of_range; + continue; + } + if (p[x] == 0.0) + ++r->zero_hits; + ++r->counts[x]; + } + free(cdf); + return r->out_of_range == 0 && sum_counts(r->counts, n) == samples; +} + +static struct gof_result goodness_of_fit(const double *p, + const uint64_t *obs, + uint32_t n, + uint64_t samples) +{ + struct gof_result r; + double tail_e = 0.0; + double tail_o = 0.0; + uint32_t bins = 0; + uint32_t i; + + r.pearson = 0.0; + r.g = 0.0; + for (i = 0; i < n; ++i) { + double e = p[i] * (double)samples; + double o = (double)obs[i]; + + if (e >= MIN_EXPECTED) { + double d = o - e; + r.pearson += d * d / e; + if (o > 0.0) + r.g += 2.0 * o * log(o / e); + ++bins; + } + else { + tail_e += e; + tail_o += o; + } + } + /* + * Do not discard the pooled small-expected tail. Dropping it breaks + * the equality between total observed and total expected counts over + * the retained bins, and can make the likelihood-ratio statistic + * negative, which then produces a NaN p-value. + */ + if (tail_e > 0.0 || tail_o > 0.0) { + if (tail_e > 0.0) { + double d = tail_o - tail_e; + r.pearson += d * d / tail_e; + if (tail_o > 0.0) + r.g += 2.0 * tail_o * log(tail_o / tail_e); + ++bins; + } + else { + r.pearson = INFINITY; + r.g = INFINITY; + ++bins; + } + } + if (r.pearson < 0.0 && r.pearson > -1.0e-10) + r.pearson = 0.0; + if (r.g < 0.0 && r.g > -1.0e-10) + r.g = 0.0; + + r.df = bins > 1 ? (int)bins - 1 : 0; + if (r.df == 0) { + r.pearson_p = 1.0; + r.g_p = 1.0; + } + else { + r.pearson_p = chi_square_pvalue(r.pearson, r.df); + r.g_p = chi_square_pvalue(r.g, r.df); + } + return r; +} + +static struct two_sample_result two_sample_chi_square(const uint64_t *a, + uint64_t na, + const uint64_t *b, + uint64_t nb, + uint32_t n) +{ + struct two_sample_result r; + double tail_a = 0.0; + double tail_b = 0.0; + uint32_t bins = 0; + uint32_t i; + + r.chi2 = 0.0; + r.max_diff = 0.0; + r.df = 0; + if (na == 0 || nb == 0) { + r.pvalue = NAN; + return r; + } + for (i = 0; i < n; ++i) { + double ai = (double)a[i]; + double bi = (double)b[i]; + double total = ai + bi; + double diff = fabs(ai / (double)na - bi / (double)nb); + + if (diff > r.max_diff) + r.max_diff = diff; + if (total > 0.0) { + double ea = (double)na * total / (double)(na + nb); + double eb = (double)nb * total / (double)(na + nb); + + if (ea >= MIN_EXPECTED && eb >= MIN_EXPECTED) { + double da = ai - ea; + double db = bi - eb; + + r.chi2 += da * da / ea + db * db / eb; + ++bins; + } + else { + tail_a += ai; + tail_b += bi; + } + } + } + if (tail_a + tail_b > 0.0) { + double total = tail_a + tail_b; + double ea = (double)na * total / (double)(na + nb); + double eb = (double)nb * total / (double)(na + nb); + + if (ea >= MIN_EXPECTED && eb >= MIN_EXPECTED) { + double da = tail_a - ea; + double db = tail_b - eb; + + r.chi2 += da * da / ea + db * db / eb; + ++bins; + } + } + if (r.chi2 < 0.0 && r.chi2 > -1.0e-10) + r.chi2 = 0.0; + + r.df = bins > 1 ? (int)bins - 1 : 0; + if (r.df == 0) + r.pvalue = 1.0; + else + r.pvalue = chi_square_pvalue(r.chi2, r.df); + return r; +} + +static int support_ok(const struct sample_result *r) +{ + return r->out_of_range == 0 && r->zero_hits == 0; +} + +static int make_uniform(double *w, uint32_t n) +{ + uint32_t i; + for (i = 0; i < n; ++i) + w[i] = 1.0; + return 1; +} + +static int make_linear(double *w, uint32_t n) +{ + uint32_t i; + for (i = 0; i < n; ++i) + w[i] = (double)(i + 1); + return 1; +} + +static int make_alternating_zeros(double *w, uint32_t n) +{ + uint32_t i; + for (i = 0; i < n; ++i) + w[i] = (i & 1U) ? 0.0 : 1.0; + return 1; +} + +static int make_dynamic_range(double *w, uint32_t n) +{ + uint32_t i; + for (i = 0; i < n; ++i) + w[i] = exp(-0.75 * (double)i); + return 1; +} + +static int make_random_weights(rng_state *rng, double *w, uint32_t n) +{ + uint32_t i; + int positive = 0; + + for (i = 0; i < n; ++i) { + uint32_t r = mt19937_next32(rng); + double u = u32_to_open01(r); + double z = u32_to_open01(mt19937_next32(rng)); + + if ((r & UINT32_C(15)) == 0) { + w[i] = 0.0; + } + else { + w[i] = exp(-30.0 * u) * (1.0 + 1000.0 * z); + positive = 1; + } + } + if (!positive) + w[0] = 1.0; + return 1; +} + +static int deterministic_distribution(const double *p, uint32_t n, + uint32_t *value) +{ + uint32_t i; + uint32_t seen = 0; + + *value = UINT32_MAX; + for (i = 0; i < n; ++i) { + if (p[i] > 0.0) { + *value = i; + ++seen; + } + } + return seen == 1; +} + +static int run_distribution_case(const char *name, const double *w, + uint32_t n, uint64_t samples, + uint32_t seed_base, + struct fisher_accumulator *fisher, + int random_case) +{ + double *p; + rng_state *vose_rng = NULL; + rng_state *wheel_rng = NULL; + rng_alias_sampler_t *sampler = NULL; + struct sample_result vose; + struct sample_result wheel; + struct gof_result gv; + struct gof_result gw; + struct two_sample_result diff; + uint32_t det_value; + int ok = 1; + int allocated_vose = 0; + int allocated_wheel = 0; + double cutoff = random_case ? RANDOM_INDIVIDUAL_MIN : STAT_P_MIN; + + p = (double *)malloc((size_t)n * sizeof(*p)); + if (p == NULL) + return 0; + if (!normalize_weights(w, p, n)) { + free(p); + return 0; + } + + sampler = rng_alias_sampler_new(w, n); + vose_rng = rng_new(seed_base + UINT32_C(1)); + wheel_rng = rng_new(seed_base + UINT32_C(2)); + if (sampler == NULL || vose_rng == NULL || wheel_rng == NULL) { + ok = 0; + goto cleanup; + } + if (!alloc_sample_result(&vose, n) || !alloc_sample_result(&wheel, n)) { + ok = 0; + goto cleanup; + } + allocated_vose = 1; + allocated_wheel = 1; + + if (!sample_vose(vose_rng, sampler, p, n, samples, &vose)) + ok = 0; + if (!sample_wheel(wheel_rng, w, p, n, samples, &wheel)) + ok = 0; + if (!support_ok(&vose) || !support_ok(&wheel)) + ok = 0; + + if (deterministic_distribution(p, n, &det_value)) { + if (vose.counts[det_value] != samples) + ok = 0; + if (wheel.counts[det_value] != samples) + ok = 0; + if (!ok) { + printf(" %-30s deterministic support FAIL\n", name); + printf(" %-30s Vose count=%llu Wheel count=%llu\n", + "", + (unsigned long long)vose.counts[det_value], + (unsigned long long)wheel.counts[det_value]); + } + goto cleanup; + } + + gv = goodness_of_fit(p, vose.counts, n, samples); + gw = goodness_of_fit(p, wheel.counts, n, samples); + diff = two_sample_chi_square(wheel.counts, samples, vose.counts, + samples, n); + + if (!isfinite(gv.pearson_p) || !isfinite(gv.g_p) || + gv.pearson_p <= cutoff || gv.g_p <= cutoff) + ok = 0; + if (!isfinite(gw.pearson_p) || gw.pearson_p <= cutoff) + ok = 0; + if (!isfinite(diff.pvalue) || diff.pvalue <= cutoff) + ok = 0; + + if (fisher != NULL) { + fisher_add(fisher, gv.pearson_p); + fisher_add(fisher, gv.g_p); + fisher_add(fisher, gw.pearson_p); + fisher_add(fisher, diff.pvalue); + } + + if (!ok) { + printf(" %-30s Vose P=%8.5f G=%8.5f Wheel P=%8.5f\n", + name, gv.pearson_p, gv.g_p, gw.pearson_p); + printf(" %-30s Diff P=%8.5f max|dp|=%8.6f FAIL\n", + "", diff.pvalue, diff.max_diff); + printf(" %-30s out_of_range V=%llu W=%llu zero_hits V=%llu W=%llu\n", + "", + (unsigned long long)vose.out_of_range, + (unsigned long long)wheel.out_of_range, + (unsigned long long)vose.zero_hits, + (unsigned long long)wheel.zero_hits); + } + +cleanup: + if (allocated_vose) + free_sample_result(&vose); + if (allocated_wheel) + free_sample_result(&wheel); + rng_alias_sampler_free(sampler); + rng_free(vose_rng); + rng_free(wheel_rng); + free(p); + return ok; +} + +static int serial_independence_case(const char *name, const double *w, + uint32_t n, uint32_t seed) +{ + double *p = NULL; + uint64_t *pairs = NULL; + rng_alias_sampler_t *sampler = NULL; + rng_state *rng = NULL; + uint32_t prev; + uint64_t i; + double stat = 0.0; + double tail_e = 0.0; + double tail_o = 0.0; + uint32_t bins = 0; + int df; + double pvalue; + int ok = 1; + + if (n > 32) + return 1; + p = (double *)malloc((size_t)n * sizeof(*p)); + pairs = (uint64_t *)calloc((size_t)n * (size_t)n, sizeof(*pairs)); + sampler = rng_alias_sampler_new(w, n); + rng = rng_new(seed); + if (p == NULL || pairs == NULL || sampler == NULL || rng == NULL || + !normalize_weights(w, p, n)) { + ok = 0; + goto cleanup; + } + + prev = rng_alias_sampler_choose(rng, sampler); + if (prev >= n || p[prev] == 0.0) { + ok = 0; + goto cleanup; + } + for (i = 1; i < SERIAL_SAMPLES; ++i) { + uint32_t cur = rng_alias_sampler_choose(rng, sampler); + if (cur >= n || p[cur] == 0.0) { + ok = 0; + goto cleanup; + } + ++pairs[(size_t)prev * n + cur]; + prev = cur; + } + + for (i = 0; i < (uint64_t)n * (uint64_t)n; ++i) { + uint32_t a = (uint32_t)(i / n); + uint32_t b = (uint32_t)(i % n); + double e = (double)(SERIAL_SAMPLES - 1) * p[a] * p[b]; + double o = (double)pairs[i]; + + if (e >= MIN_EXPECTED) { + double d = o - e; + stat += d * d / e; + ++bins; + } + else { + tail_e += e; + tail_o += o; + } + } + if (tail_e > 0.0 || tail_o > 0.0) { + if (tail_e > 0.0) { + double d = tail_o - tail_e; + stat += d * d / tail_e; + ++bins; + } + else { + stat = INFINITY; + ++bins; + } + } + df = bins > 1 ? (int)bins - 1 : 0; + pvalue = chi_square_pvalue(stat, df); + if (!isfinite(pvalue) || pvalue <= STAT_P_MIN) + ok = 0; + + printf(" %-30s serial P=%8.5f %s\n", + name, pvalue, ok ? "PASS" : "FAIL"); + +cleanup: + free(pairs); + free(p); + rng_alias_sampler_free(sampler); + rng_free(rng); + return ok; +} + +static int test_invalid_inputs(void) +{ + rng_alias_sampler_t *s; + double p2[] = { 0.5, 0.5 }; + double zero[] = { 0.0, 0.0 }; + double negative[] = { 0.5, -0.1, 0.6 }; + double nanv[] = { 0.5, NAN }; + int ok = 1; + + s = rng_alias_sampler_new(NULL, 0); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + s = rng_alias_sampler_new(NULL, 2); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + s = rng_alias_sampler_new(p2, 0); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + s = rng_alias_sampler_new(zero, 2); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + s = rng_alias_sampler_new(negative, 3); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + s = rng_alias_sampler_new(nanv, 2); + if (s != NULL) { rng_alias_sampler_free(s); ok = 0; } + return ok; +} + +static int test_lifecycle(void) +{ + double w[] = { 7.0, 0.0, 3.0, 1.0, 11.0 }; + uint32_t i; + + for (i = 0; i < LIFECYCLE_ITERATIONS; ++i) { + rng_alias_sampler_t *s = rng_alias_sampler_new(w, 5); + if (s == NULL) + return 0; + rng_alias_sampler_free(s); + } + return 1; +} + +static int run_named_distribution_tests(void) +{ + struct fisher_accumulator fisher; + int ok = 1; + double *w; + + fisher_init(&fisher); + w = (double *)malloc((size_t)4099 * sizeof(*w)); + if (w == NULL) + return 0; + +#define RUN_CASE(label, nval, maker, samples, seed) \ + do { \ + if (!(maker)(w, (nval)) || \ + !run_distribution_case((label), w, (nval), (samples), \ + (seed), &fisher, 0)) \ + ok = 0; \ + } while (0) + + RUN_CASE("n=1", 1, make_uniform, STANDARD_SAMPLES, 1000); + RUN_CASE("n=2 half", 2, make_uniform, STANDARD_SAMPLES, 2000); + RUN_CASE("uniform n=3", 3, make_uniform, STANDARD_SAMPLES, 3000); + RUN_CASE("uniform n=7", 7, make_uniform, STANDARD_SAMPLES, 4000); + RUN_CASE("linear n=17", 17, make_linear, STANDARD_SAMPLES, 5000); + RUN_CASE("zeros n=31", 31, make_alternating_zeros, + STANDARD_SAMPLES, 6000); + RUN_CASE("dynamic n=64", 64, make_dynamic_range, + STANDARD_SAMPLES, 7000); + RUN_CASE("large n=4099", 4099, make_uniform, + LARGE_TABLE_SAMPLES, 8000); + +#undef RUN_CASE + + { + double boundary2[] = { 0.5 - DBL_EPSILON, 0.5 + DBL_EPSILON }; + double near_zero[] = { 1.0, DBL_EPSILON, 0.0, DBL_MIN }; + double near_one[] = { 1.0 - 1.0e-9, 1.0e-9, 0.0, 0.0 }; + double unnorm[] = { 1000.0, 0.0, 1.0, 3.0, 11.0, 31.0 }; + double extreme[] = { 0.0, 1e-300, 1e-200, 1e-100, + 1e-50, 1e-20, 1e-10, 1.0, + 100.0, 1e5, 1e10 }; + + if (!run_distribution_case("boundary 0.5", boundary2, 2, + STANDARD_SAMPLES, 9000, &fisher, 0)) + ok = 0; + if (!run_distribution_case("near zero", near_zero, 4, + STANDARD_SAMPLES, 10000, &fisher, 0)) + ok = 0; + if (!run_distribution_case("near one", near_one, 4, + STANDARD_SAMPLES, 11000, &fisher, 0)) + ok = 0; + if (!run_distribution_case("unnormalized", unnorm, 6, + STANDARD_SAMPLES, 12000, &fisher, 0)) + ok = 0; + if (!run_distribution_case("extreme range", extreme, 11, + STANDARD_SAMPLES, 13000, &fisher, 0)) + ok = 0; + if (!serial_independence_case("serial uniform n=7", w, 7, + MASTER_SEED + 14000)) + ok = 0; + if (!serial_independence_case("serial boundary 0.5", boundary2, 2, + MASTER_SEED + 15000)) + ok = 0; + } + { + double fp = fisher_pvalue(&fisher); + printf(" %-30s Fisher P=%8.5f (%u tests)\n", + "standard aggregate", fp, fisher.count); + if (!isfinite(fp) || fp <= FISHER_P_MIN) + ok = 0; + } + free(w); + return ok; +} + +static int run_randomized_tests(void) +{ + struct fisher_accumulator fisher; + rng_state *master; + double *w; + uint32_t case_no; + int ok = 1; + + fisher_init(&fisher); + master = rng_new(MASTER_SEED + UINT32_C(0x12345678)); + w = (double *)malloc((size_t)RANDOM_MAX_N * sizeof(*w)); + if (master == NULL || w == NULL) { + rng_free(master); + free(w); + return 0; + } + for (case_no = 0; case_no < RANDOM_CASES; ++case_no) { + uint32_t n = 2 + rng_uniform_u32_ubound(master, RANDOM_MAX_N - 1); + uint32_t seed = MASTER_SEED + UINT32_C(200000) + case_no * 4U; + int case_ok; + + if (!make_random_weights(master, w, n)) { + ok = 0; + break; + } + case_ok = run_distribution_case("random", w, n, RANDOM_SAMPLES, + seed, &fisher, 1); + if (!case_ok) + ok = 0; + } + { + double fp = fisher_pvalue(&fisher); + printf(" %-30s Fisher P=%8.5f (%u tests)\n", + "random aggregate", fp, fisher.count); + if (!isfinite(fp) || fp <= FISHER_P_MIN) + ok = 0; + } + rng_free(master); + free(w); + return ok; +} + +int main(void) +{ + int ok = 1; + int result; + + printf("Walker-Vose alias sampler test harness\n"); + printf("======================================\n\n"); + + result = test_invalid_inputs(); + record_result("invalid inputs", result); + if (!result) + ok = 0; + + result = test_lifecycle(); + record_result("repeated construction/destruction", result); + if (!result) + ok = 0; + + printf("\nNamed and boundary distributions\n"); + printf("--------------------------------\n"); + result = run_named_distribution_tests(); + record_result("named distribution suite", result); + if (!result) + ok = 0; + + printf("\nRandomized distributions\n"); + printf("------------------------\n"); + result = run_randomized_tests(); + record_result("randomized distribution suite", result); + if (!result) + ok = 0; + + printf("\n======================================\n"); + if (ok) { + printf("RESULT: PASS (%u test groups)\n", g_tests_run); + } + else { + printf("RESULT: FAIL (%u test groups, %u failed)\n", + g_tests_run, g_tests_failed); + } + return ok ? EXIT_SUCCESS : EXIT_FAILURE; +}