FPRAS

.github/workflows/R.yml

@@ -190,8 +190,9 @@ jobs:
     needs: R-CMD-check
     permissions:
       contents: write # github-pages-deploy-action
+      pull-requests: write # post PR preview comment
     concurrency: # Recommended if you intend to make multiple deployments in quick succession.
-      group: web-${{ github.workflow }}-${{ github.ref }}
+      group: web-${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -233,6 +234,46 @@ jobs:
           single-commit: true
           clean: true
 
+      - name: Deploy PR preview 🚀
+        if: success() && github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository
+        uses: JamesIves/github-pages-deploy-action@v4
+        with:
+          folder: r/docs
+          target-folder: previews/pr-${{ github.event.pull_request.number }}/r
+          single-commit: true
+          clean: true
+
+      - name: Comment PR preview link
+        if: success() && github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository
+        uses: actions/github-script@v7
+        with:
+          script: |
+            const marker = '<!-- docs-preview-r-marker -->';
+            const prNumber = context.issue.number;
+            const url = `https://multi-objective.github.io/moocore/previews/pr-${prNumber}/r/`;
+            const body = `${marker}\n📚 [**R docs preview** for this PR](${url})\n\n`;
+            const comments = await github.rest.issues.listComments({
+              owner: context.repo.owner,
+              repo: context.repo.repo,
+              issue_number: prNumber,
+            });
+            const existing = comments.data.find(c => c.body.includes(marker));
+            if (existing) {
+              await github.rest.issues.updateComment({
+                owner: context.repo.owner,
+                repo: context.repo.repo,
+                comment_id: existing.id,
+                body,
+              });
+            } else {
+              await github.rest.issues.createComment({
+                owner: context.repo.owner,
+                repo: context.repo.repo,
+                issue_number: prNumber,
+                body,
+              });
+            }
+
   rchk:
     needs: pkgdown
     uses: multi-objective/moocore/.github/workflows/rchk.yml@main

.github/workflows/docs-preview-cleanup.yml

@@ -0,0 +1,34 @@
+name: Cleanup docs preview
+
+on:
+  pull_request:
+    types: [closed]
+
+permissions:
+  contents: write
+
+jobs:
+  cleanup:
+    name: Remove docs preview for closed PR
+    runs-on: ubuntu-latest
+    # Only run for same-repo PRs; fork PRs never get a preview deployed
+    if: github.event.pull_request.head.repo.full_name == github.repository
+    steps:
+      - uses: actions/checkout@v6
+        with:
+          ref: gh-pages
+
+      - name: Remove preview folder
+        run: |
+          PR_NUMBER="${{ github.event.pull_request.number }}"
+          PREVIEW_DIR="previews/pr-${PR_NUMBER}"
+          if [ -d "${PREVIEW_DIR}" ]; then
+            git rm -rf "${PREVIEW_DIR}"
+            git config user.email "actions@github.com"
+            git config user.name "GitHub Actions"
+            git commit -m "Remove docs preview for PR #${PR_NUMBER}"
+            git push
+            echo "Removed ${PREVIEW_DIR}"
+          else
+            echo "No preview directory found for PR #${PR_NUMBER}, nothing to clean up"
+          fi

.github/workflows/python.yml

@@ -117,8 +117,9 @@ jobs:
     if: "! contains(github.event.head_commit.message, '[skip ci]')"
     permissions:
       contents: write
+      pull-requests: write # post PR preview comment
     concurrency: # Recommended if you intend to make multiple deployments in quick succession.
-      group: web-${{ github.workflow }}-${{ github.ref }}
+      group: web-${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
     timeout-minutes: 10
     runs-on:  ubuntu-latest
     steps:
@@ -165,6 +166,46 @@ jobs:
           single-commit: true
           clean: true
 
+      - name: Deploy PR preview 🚀
+        if: success() && github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository
+        uses: JamesIves/github-pages-deploy-action@v4
+        with:
+          folder: python/doc/_build/html/
+          target-folder: previews/pr-${{ github.event.pull_request.number }}/python
+          single-commit: true
+          clean: true
+
+      - name: Comment PR preview link
+        if: success() && github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository
+        uses: actions/github-script@v7
+        with:
+          script: |
+            const marker = '<!-- docs-preview-python-marker -->';
+            const prNumber = context.issue.number;
+            const url = `https://multi-objective.github.io/moocore/previews/pr-${prNumber}/python/`;
+            const body = `${marker}\n📚 [**Python docs preview** for this PR](${url})\n\n`;
+            const comments = await github.rest.issues.listComments({
+              owner: context.repo.owner,
+              repo: context.repo.repo,
+              issue_number: prNumber,
+            });
+            const existing = comments.data.find(c => c.body.includes(marker));
+            if (existing) {
+              await github.rest.issues.updateComment({
+                owner: context.repo.owner,
+                repo: context.repo.repo,
+                comment_id: existing.id,
+                body,
+              });
+            } else {
+              await github.rest.issues.createComment({
+                owner: context.repo.owner,
+                repo: context.repo.repo,
+                issue_number: prNumber,
+                body,
+              });
+            }
+
   coverage:
     timeout-minutes: 15
     needs: build-doc

c/config.h

@@ -25,6 +25,13 @@ typedef uint_fast8_t dimension_t;
 // Maximum number of objectives is 31. It cannot be larger than 255.
 #define MOOCORE_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/hvapprox.c

@@ -5,6 +5,7 @@
 #include "hvapprox.h"
 #include "pow_int.h"
 #include "rng.h"
+#include "sort.h"
 
 #define ALMOST_ZERO_WEIGHT 1e-20
 
@@ -863,3 +864,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;
+}
+
+static double
+cumsum_of_vector(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];
+}
+
+/**
+   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.
+
+   delta : Probability of failure.
+*/
+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_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)
+    // Convert vols into a piece-wise cdf.
+    double total_vol = cumsum_of_vector(vols, npoints);
+    for (size_t k = 0; k < npoints; k++)
+        vols[k] /= total_vol;
+
+    rng_state * rng = rng_new(random_seed);
+    uint64_t t_sum = 0, m  = 0;
+    while (true) {
+        // FIXME: Ideally, this should be done in O(1), but currently it takes O(log n).
+        uint32_t i = rng_random_wheel_uint32(rng, vols, (uint32_t) npoints);
+
+        // 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);
+                free((void *) points);
+                return total_vol * (T_factor / (double) m);
+            }
+            t_sum++;
+            uint32_t j = rng_uint32(rng, 0, (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

@@ -34,6 +34,49 @@ rng_uniform(rng_state * rng, double low, double high)
     return low + (high - low) * rng_random(rng);
 }
 
+static inline uint32_t
+rng_uint32(rng_state * rng, uint32_t low, uint32_t high)
+{
+    assert(rng != NULL);
+    if (low >= high)
+        return low;
+    return low + (mt19937_next32(rng) % (high - low));
+}
+
+/**
+   Returns a random value within [0, n - 1] based on the given CDF.
+
+   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)
+{
+    // Check the CDF is correct.
+    DEBUG1(for (uint32_t j = 1; j < n; j++) assert(cdf[j-1] <= cdf[j]));
+    DEBUG1(for (uint32_t j = 0; j < n; j++) assert(0 <= cdf[j] && cdf[j] <= 1));
+
+    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,

python/doc/source/reference/functions.metrics.rst

@@ -201,6 +201,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