Add hv_approx_fpras() in R (FPRAS hypervolume approximation)

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

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_ */

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

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 */

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.
+
+   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_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);
+                if (unlikely(m == 0)) return 0.0;
+                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++;
+    }
+}

c/hvapprox.h

@@ -37,5 +37,17 @@ 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.
+*/
+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_

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/32306366
+
+*/
+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 */