FPRAS

bibkeys.txt

@@ -41,6 +41,7 @@ RubMel1998simulation
 SchEsqLarCoe2012tec
 SchKer2025r2v2
 VelLam1998gp
+Vose1991alias
 WuAza2001metrics
 ZhoZhaJin2009igdx
 ZitThi1998ppsn

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.
+   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++;
+    }
+}

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_