diff options
| -rw-r--r-- | CMakeLists.txt | 20 | ||||
| -rw-r--r-- | src/extent-sampling.c | 161 | ||||
| -rw-r--r-- | src/extent-sampling.sc | 155 | ||||
| -rw-r--r-- | src/gaussian-nd-random.c | 51 | ||||
| -rw-r--r-- | src/gaussian-nd-random.sc | 72 | ||||
| -rw-r--r-- | src/macros/macros.sc | 12 | ||||
| -rw-r--r-- | src/nd-random.c | 101 | ||||
| -rw-r--r-- | src/nd-random.sc | 98 | ||||
| -rw-r--r-- | src/oracles.c | 108 | ||||
| -rw-r--r-- | src/oracles.sc | 90 | ||||
| -rw-r--r-- | src/utils.c | 160 | ||||
| -rw-r--r-- | src/utils.sc | 123 |
12 files changed, 568 insertions, 583 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 94264c7..fc3a296 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,8 +2,24 @@ cmake_minimum_required(VERSION 3.10) project(extent-sampling VERSION 0.1) find_package(GSL REQUIRED) +# TODO: Make indent optional +find_program(INDENT NAMES indent REQUIRED) +find_program(SC NAMES sc REQUIRED) + +# Generate C source files from SC source files. +file(GLOB SC_SOURCES "src/*.sc") +foreach(sc_source ${SC_SOURCES}) + get_filename_component(source_basename ${sc_source} NAME_WE) + set(c_file ${source_basename}.c) + add_custom_command( + OUTPUT ${c_file} + COMMAND ${SC} ${sc_source} ${c_file} + COMMAND ${INDENT} ${c_file} + DEPENDS ${sc_source} + VERBATIM) + list(APPEND C_SOURCES ${c_file}) +endforeach() include_directories("include") -file(GLOB SOURCES "src/*.c") -add_library(extentsampling SHARED ${SOURCES}) +add_library(extentsampling SHARED ${C_SOURCES}) diff --git a/src/extent-sampling.c b/src/extent-sampling.c deleted file mode 100644 index eff3b3f..0000000 --- a/src/extent-sampling.c +++ /dev/null @@ -1,161 +0,0 @@ -/* - This file contains all the core extent sampling routines. From time - to time as research progresses, I keep refactoring this code to be - clearer and more efficient. This mostly involves throwing out old - code even when such old code may be required in the future. I - figured that keeping old code around is not worth the trouble - because it impedes clarity of thought. -*/ - -#include <math.h> -#include <stddef.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_integration.h> -#include <gsl/gsl_rstat.h> -#include <gsl/gsl_randist.h> -#include "extent-sampling.h" -#include "nd-random.h" -#include "utils.h" - -#define VOLUME_MINIMUM_SAMPLES 100 -#define INTEGRAL_MINIMUM_SAMPLES 100 -#define CONFIDENCE_INTERVAL_FACTOR 1.96 - -#define INTEGRATION_INTERVALS 1000 - -#define SOLID_ANGLE_START 0.0 -#define SOLID_ANGLE_LAST 0.5 -#define SOLID_ANGLE_STEPS 100 - -double volume (extent_oracle_t extent_oracle, double true_volume, - const gsl_rng* r, unsigned int dimension, double rtol, - gsl_rstat_workspace* stats) -{ - gsl_vector* x = gsl_vector_alloc(dimension); - double volume, error; - double vn = ln_volume_of_ball(dimension); - do { - random_direction_vector(r, x); - double extent = extent_oracle(x); - double volume_from_sample = exp(vn + dimension*log(extent)); - gsl_rstat_add(volume_from_sample, stats); - volume = gsl_rstat_mean(stats); - error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / volume; - } while ((error > rtol) || (rerror(volume, true_volume) > rtol) || (gsl_rstat_n(stats) < VOLUME_MINIMUM_SAMPLES)); - gsl_vector_free(x); - return volume; -} - -double volume_window (extent_oracle_t extent_oracle, double true_volume, - const gsl_rng* r, unsigned int dimension, double rtol, - unsigned int* number_of_samples) -{ - gsl_rstat_workspace* stats = gsl_rstat_alloc(); - gsl_vector* x = gsl_vector_alloc(dimension); - double volume, error; - double vn = ln_volume_of_ball(dimension); - // This is the window length used in Volume.m of Lovasz-Vempala's - // code. - int window_length = 4*dimension*dimension + 500; - int accurate_estimates = 0; - do { - random_direction_vector(r, x); - double extent = extent_oracle(x); - double volume_from_sample = exp(vn + dimension*log(extent)); - gsl_rstat_add(volume_from_sample, stats); - volume = gsl_rstat_mean(stats); - error = rerror(volume, true_volume); - if (error < rtol) accurate_estimates++; - else accurate_estimates = 0; - } while (accurate_estimates < window_length); - *number_of_samples = gsl_rstat_n(stats); - gsl_vector_free(x); - gsl_rstat_free(stats); - return volume; -} - -double integral_per_direction -(integrand_t integrand, const gsl_vector* direction, const gsl_rng* r, unsigned int n, - double radius, double rtol, int* neval) -{ - double f (double r, void* params) - { - (*neval)++; - return gsl_pow_uint(r, n-1) * integrand(r, direction); - } - - gsl_function gsl_f = {&f}; - double result, error; - gsl_integration_workspace *w = gsl_integration_workspace_alloc(INTEGRATION_INTERVALS); - gsl_integration_qag(&gsl_f, 0, radius, 0, rtol, INTEGRATION_INTERVALS, GSL_INTEG_GAUSS15, w, &result, &error); - gsl_integration_workspace_free(w); - return surface_area_of_ball(n) * result; -} - -double integral -(integrand_t integrand, extent_oracle_t extent_oracle, double true_integral, - const gsl_rng* r, unsigned int dimension, double rtol, - gsl_rstat_workspace* stats) -{ - gsl_vector* x = gsl_vector_alloc(dimension); - double integral, error; - do { - int neval = 0; - random_direction_vector(r, x); - double extent = extent_oracle(x); - double integral_from_sample = integral_per_direction(integrand, x, r, dimension, extent, rtol, &neval); - gsl_rstat_add(integral_from_sample, stats); - integral = gsl_rstat_mean(stats); - error = CONFIDENCE_INTERVAL_FACTOR * gsl_rstat_sd_mean(stats) / integral; - } while ((error > rtol) || (rerror(integral, true_integral) > rtol) || (gsl_rstat_n(stats) < INTEGRAL_MINIMUM_SAMPLES)); - if (error > rtol) - GSL_ERROR_VAL("integral failed to converge", GSL_ETOL, integral); - gsl_vector_free(x); - return integral; -} - -double volume_cone (extent_oracle_t extent_oracle, const gsl_rng* r, - const gsl_vector* mean, double omega_min, double omega_max, - unsigned int number_of_samples, double* variance) -{ - unsigned int n = mean->size; - gsl_vector* x = gsl_vector_alloc(n); - double theta_min = solid_angle_to_planar_angle(omega_min * surface_area_of_ball(n), n); - double theta_max = solid_angle_to_planar_angle(omega_max * surface_area_of_ball(n), n); - double vn = ln_volume_of_ball(n); - gsl_rstat_workspace* stats = gsl_rstat_alloc(); - for (int i=0; i<number_of_samples; i++) { - hollow_cone_random_vector(r, mean, theta_min, theta_max, x); - double extent = extent_oracle(x); - double volume_from_sample = exp(vn + n*log(extent)); - gsl_rstat_add(volume_from_sample, stats); - } - - if (variance) *variance = gsl_rstat_variance(stats); - - gsl_rstat_free(stats); - gsl_vector_free(x); - return gsl_rstat_mean(stats) * (omega_max - omega_min); -} - -double volume_experiment -(extent_oracle_t extent_oracle, const gsl_rng* r, - const gsl_vector* mean, unsigned int samples_per_cone, - double solid_angle_factor, double solid_angle_threshold_exponent_factor, - unsigned int* number_of_samples) -{ - int n = mean->size; - double volume = 0; - double omega_max = 0.5; - int N = 0; - do { - double omega_min = omega_max / solid_angle_factor; - volume += volume_cone(extent_oracle, r, mean, omega_min, omega_max, samples_per_cone, NULL); - N += samples_per_cone; - omega_max = omega_min; - } while (omega_max > pow(2, -solid_angle_threshold_exponent_factor*n)); - volume += volume_cone(extent_oracle, r, mean, 0, omega_max, samples_per_cone, NULL); - N += samples_per_cone; - if (number_of_samples) *number_of_samples = N; - return volume; -} diff --git a/src/extent-sampling.sc b/src/extent-sampling.sc new file mode 100644 index 0000000..c6220e6 --- /dev/null +++ b/src/extent-sampling.sc @@ -0,0 +1,155 @@ +(sc-include "macros/macros") + +(pre-include "stddef.h") +(pre-include "gsl/gsl_integration.h") +(pre-include "gsl/gsl_rstat.h") +(pre-include "gsl/gsl_vector.h") +(pre-include "nd-random.h") +(pre-include "utils.h") +(pre-include "extent-sampling.h") + +(sc-define-syntax (with-vector var size body ...) + (with-alloc var gsl-vector* + (gsl-vector-alloc size) + gsl-vector-free + body ...)) + +(sc-define-syntax (with-rstats var body ...) + (with-alloc var gsl-rstat-workspace* + (gsl-rstat-alloc) gsl-rstat-free + body ...)) + +(sc-define-syntax (with-integration-workspace var intervals body ...) + (with-alloc var gsl-integration-workspace* + (gsl-integration-workspace-alloc intervals) + gsl-integration-workspace-free + body ...)) + +(pre-define CONFIDENCE-INTERVAL-FACTOR 1.96) + +(pre-let* (VOLUME-MINIMUM-SAMPLES 100) + (define (volume extent-oracle true-volume r dimension rtol stats) + (double extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*) + (declare volume double) + (let* ((vn double (ln-volume-of-ball dimension))) + (with-vector x dimension + (do-while (or (> (/ (* CONFIDENCE-INTERVAL-FACTOR + (gsl-rstat-sd-mean stats)) + volume) + rtol) + (> (rerror volume true-volume) rtol) + (< (gsl-rstat-n stats) VOLUME-MINIMUM-SAMPLES)) + (random-direction-vector r x) + (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x))))) + stats) + (set volume (gsl-rstat-mean stats))))) + (return volume))) + +(define (volume-window extent-oracle true-volume r dimension rtol number-of-samples) + (double extent-oracle-t double (const gsl-rng*) (unsigned int) double (unsigned int*)) + (declare volume double) + (let* ((vn double (ln-volume-of-ball dimension)) + ;; This is the window length used in Volume.m of + ;; Lovasz-Vempala's code. + (window-length int (+ (* 4 dimension dimension) 500)) + (accurate-estimates int 0)) + (with-rstats stats + (with-vector x dimension + (do-while (< accurate-estimates window-length) + (random-direction-vector r x) + (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x))))) + stats) + (set volume (gsl-rstat-mean stats)) + (cond + ((< (rerror volume true-volume) rtol) + (set+ accurate-estimates 1)) + (else (set accurate-estimates 0))))) + (set (pointer-get number-of-samples) + (gsl-rstat-n stats)))) + (return volume)) + +(pre-let* (INTEGRATION-INTERVALS 1000) + (define (integral-per-direction integrand direction r n radius rtol neval) + (double integrand-t (const gsl-vector*) (const gsl-rng*) (unsigned int) double double int*) + (define (f r params) (double double void*) + (set+ (pointer-get neval) 1) + (return (* (gsl-pow-uint r (- n 1)) + (integrand r direction)))) + (declare gsl-f (struct-variable gsl-function (function (address-of f))) + result double + error double) + (with-integration-workspace w INTEGRATION-INTERVALS + (gsl-integration-qag (address-of gsl-f) + 0 + radius + 0 + rtol + INTEGRATION-INTERVALS + GSL-INTEG-GAUSS15 + w + (address-of result) + (address-of error))) + (return (* (surface-area-of-ball n) + result)))) + +(pre-let* (INTEGRAL-MINIMUM-SAMPLES 100) + (define (integral integrand extent-oracle true-integral r dimension rtol stats) + (double integrand-t extent-oracle-t double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*) + (declare integral double + error double) + (with-vector x dimension + (do-while (or (> error rtol) + (> (rerror integral true-integral) rtol) + (< (gsl-rstat-n stats) INTEGRAL-MINIMUM-SAMPLES)) + ;; TODO: If this neval is not used anywhere, remove it. + (let* ((neval int 0)) + (random-direction-vector r x) + (gsl-rstat-add (integral-per-direction integrand x r dimension + (extent-oracle x) rtol (address-of neval)) + stats)) + (set integral (gsl-rstat-mean stats)) + (set error (/ (* CONFIDENCE-INTERVAL-FACTOR (gsl-rstat-sd-mean stats)) + integral)))) + (cond + ((> error rtol) + (GSL-ERROR-VAL "Integral failed to converge" GSL-ETOL integral))) + (return integral))) + +(define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance) + (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*) + (declare volume double) + (let* ((dimension (unsigned int) (: mean size)) + (vn double (ln-volume-of-ball dimension)) + (theta-min double (solid-angle->planar-angle (* omega-min (surface-area-of-ball dimension)) dimension)) + (theta-max double (solid-angle->planar-angle (* omega-max (surface-area-of-ball dimension)) dimension))) + (with-rstats stats + (with-vector x dimension + (for-i i number-of-samples + (hollow-cone-random-vector r mean theta-min theta-max x) + (gsl-rstat-add (exp (+ vn (* dimension (log (extent-oracle x))))) + stats))) + (cond + (variance (set (pointer-get variance) (gsl-rstat-variance stats)))) + (set volume (* (gsl-rstat-mean stats) + (- omega-max omega-min))))) + (return volume)) + +(define (volume-experiment extent-oracle r mean samples-per-cone solid-angle-factor + solid-angle-threshold-exponent-factor number-of-samples) + (double extent-oracle-t (const gsl-rng*) (const gsl-vector*) + (unsigned int) double double (unsigned int*)) + (define + volume double 0 + omega-max double 0.5 + N int 0) + (let* ((dimension int (: mean size))) + (do-while (> omega-max (pow 2 (- (* dimension solid-angle-threshold-exponent-factor)))) + (let* ((omega-min double (/ omega-max solid-angle-factor))) + (set+ volume (volume-cone extent-oracle r mean omega-min omega-max samples-per-cone NULL)) + (set+ N samples-per-cone) + (set omega-max omega-min)))) + (set+ volume (volume-cone extent-oracle r mean 0 omega-max samples-per-cone NULL)) + (set+ N samples-per-cone) + (cond + (number-of-samples (set (pointer-get number-of-samples) N))) + (return volume)) diff --git a/src/gaussian-nd-random.c b/src/gaussian-nd-random.c deleted file mode 100644 index 16abc94..0000000 --- a/src/gaussian-nd-random.c +++ /dev/null @@ -1,51 +0,0 @@ -#include <gsl/gsl_blas.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_sf_gamma.h> -#include "nd-random.h" -#include "utils.h" - -#define INTEGRATION_EPSABS 1e-3 -#define INTEGRATION_EPSREL 1e-3 - -unsigned int shifted_gaussian_random_vector -(const gsl_rng* r, const gsl_vector* mean, double theta_max, double standard_deviation, gsl_vector* x) -{ - for (int i=0; i<x->size; i++) - gsl_vector_set(x, i, gsl_ran_gaussian(r, standard_deviation)); - gsl_vector_add(x, mean); - gsl_blas_dscal(1/gsl_blas_dnrm2(x), x); - if (angle_between_vectors(x, mean) > theta_max) - return 1 + shifted_gaussian_random_vector(r, mean, theta_max, standard_deviation, x); - else return 1; -} - -double shifted_gaussian_pdf (double theta, double mean, double theta_max, - double standard_deviation, unsigned int dimension, gsl_integration_workspace* w) -{ - // NOTE: Due to the rejection sampling, the returned result is only - // proportional to the actual PDF. - if (theta > theta_max) return 0; - - double projection = mean*cos(theta); - double integrand (double r, void* params) - { - return gsl_pow_uint(r, dimension - 1) - * gaussian_pdf((r - projection) / standard_deviation); - } - - gsl_function gsl_integrand = {&integrand}; - double result, error; - /* gsl_integration_qagiu(&gsl_integrand, 0, INTEGRATION_EPSABS, INTEGRATION_EPSREL, */ - /* w->limit, w, &result, &error); */ - double rmax = 0.5 * (projection + sqrt(gsl_pow_2(projection) + 4*(dimension-1)*gsl_pow_2(standard_deviation))); - double halfwidth = 3 * standard_deviation; - gsl_integration_qag(&gsl_integrand, - rmax - halfwidth < 0 ? 0 : rmax - halfwidth, - rmax + halfwidth, - INTEGRATION_EPSABS, INTEGRATION_EPSREL, - w->limit, GSL_INTEG_GAUSS41, w, &result, &error); - return result - * gaussian_pdf(sqrt(gsl_pow_2(mean) - gsl_pow_2(projection)) / standard_deviation) - / gsl_pow_uint(M_SQRT2*sqrt(M_PI), dimension - 2) - / gsl_pow_uint(standard_deviation, dimension); -} diff --git a/src/gaussian-nd-random.sc b/src/gaussian-nd-random.sc new file mode 100644 index 0000000..bfde566 --- /dev/null +++ b/src/gaussian-nd-random.sc @@ -0,0 +1,72 @@ +(sc-include "macros/macros") + +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_integration.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "utils.h") + +(pre-define INTEGRATION-EPSABS 1e-3) +(pre-define INTEGRATION-EPSREL 1e-3) + +(define (shifted-gaussian-random-vector r mean theta-max standard-deviation x) + ((unsigned int) (const gsl-rng*) (const gsl-vector*) double double gsl-vector*) + "Return a random vector generated by the shifted Gaussian method." + (for-i i (: x size) + (gsl-vector-set x i (gsl-ran-gaussian r standard-deviation))) + (gsl-vector-add x mean) + (gsl-blas-dscal (/ 1 (gsl-blas-dnrm2 x)) x) + (cond + ((> (angle-between-vectors x mean) theta-max) + (return (+ 1 (shifted-gaussian-random-vector r mean theta-max standard-deviation x)))) + (else (return 1)))) + +(define (shifted-gaussian-pdf theta mean theta-max standard-deviation dimension w) + (double double double double double (unsigned int) gsl-integration-workspace*) + "Return the unnormalized probability density function at THETA for a +random vector generated using the shifted Gaussian method." + (cond + ((> theta theta-max) (return 0)) + (else + (define projection double (* mean (cos theta))) + (define (integrand r params) (double double void*) + (return (* (gsl-pow-uint r (- dimension 1)) + (gaussian-pdf (/ (- r projection) + standard-deviation))))) + + (declare gsl-integrand (struct-variable gsl-function (function (address-of integrand))) + result double + error double) + + ;; (gsl-integration-qagiu (address-of gsl-integrand) + ;; 0 + ;; INTEGRATION-EPSABS + ;; INTEGRATION-EPSREL + ;; (: w limit) + ;; w + ;; (address-of result) + ;; (address-of error)) + (let* ((rmax double (/ (+ projection + (sqrt (+ (gsl-pow-2 projection) + (* 4 + (- dimension 1) + (gsl-pow-2 standard-deviation))))) + 2)) + (half-width double (* 3 standard-deviation))) + (gsl-integration-qag (address-of gsl-integrand) + (if* (< (- rmax half-width) 0) + 0 + (- rmax half-width)) + (+ rmax half-width) + INTEGRATION-EPSABS + INTEGRATION-EPSREL + (: w limit) + GSL-INTEG-GAUSS41 + w + (address-of result) + (address-of error))) + (return (/ (* result + (gaussian-pdf (/ (sqrt (- (gsl-pow-2 mean) + (gsl-pow-2 projection))) + standard-deviation))) + (gsl-pow-uint (sqrt (* 2 M-PI)) (- dimension 2)) + (gsl-pow-uint standard-deviation dimension)))))) diff --git a/src/macros/macros.sc b/src/macros/macros.sc new file mode 100644 index 0000000..6ffc7ec --- /dev/null +++ b/src/macros/macros.sc @@ -0,0 +1,12 @@ +(sc-define-syntax (for-i index limit body ...) + (for ((define index int 0) + (< index limit) + (set+ index 1)) + body ...)) + +(sc-define-syntax (with-alloc var type allocate free body ...) + (let* () + (declare var type) + (set var allocate) + body ... + (free var))) diff --git a/src/nd-random.c b/src/nd-random.c deleted file mode 100644 index 25abc7a..0000000 --- a/src/nd-random.c +++ /dev/null @@ -1,101 +0,0 @@ -#include <float.h> -#include <math.h> -#include <gsl/gsl_blas.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_sf_gamma.h> -#include "nd-random.h" -#include "utils.h" - -static double beta_inc_unnormalized (double a, double b, double x); -static double incomplete_wallis_integral (double theta, unsigned int m); - -void random_direction_vector (const gsl_rng* r, gsl_vector* x) -{ - gsl_ran_dir_nd(r, x->size, x->data); -} - -static void rotate_from_nth_canonical (gsl_vector* x, const gsl_vector* orient) -{ - const size_t n = x->size; - double xn = gsl_vector_get(x, n - 1); - double mun = gsl_vector_get(orient, n - 1); - gsl_vector_const_view orient_sub = gsl_vector_const_subvector(orient, 0, n - 1); - double b = gsl_blas_dnrm2(&orient_sub.vector); - double a = (dot_product(orient, x) - xn*mun) / b; - double c = mun, s = sqrt(1 - gsl_pow_2(c)); - gsl_blas_daxpy((xn*s + a*(c - 1))/b, orient, x); - gsl_vector_set(x, n - 1, gsl_vector_get(x, n - 1) - + xn*(c - 1) - a*s - - mun*(xn*s + a*(c - 1))/b); -} - -/* TODO: There is an edge case when mean is the (n-1)th canonical - basis vector. Fix it. -*/ -void hollow_cone_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_min, double theta_max, gsl_vector* x) -{ - unsigned int n = x->size; - gsl_ran_dir_nd(r, n - 1, x->data); - - // Generate random vector around the nth canonical basis vector - double omega_min = planar_angle_to_solid_angle(theta_min, n); - double omega_max = planar_angle_to_solid_angle(theta_max, n); - gsl_vector_scale(x, cos(theta_max) * tan(solid_angle_to_planar_angle(gsl_ran_flat(r, omega_min, omega_max), n))); - gsl_vector_set(x, n - 1, cos(theta_max)); - gsl_vector_scale(x, 1.0/gsl_blas_dnrm2(x)); - - // Rotate to arbitrary basis - rotate_from_nth_canonical(x, mean); -} - -void subsampling_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_max, gsl_vector* x) -{ - hollow_cone_random_vector(r, mean, 0, theta_max, x); -} - -static double beta_inc_unnormalized (double a, double b, double x) -{ - return gsl_sf_beta_inc(a, b, x) * gsl_sf_beta(a, b); -} - -static double incomplete_wallis_integral (double theta, unsigned int m) -{ - /** - @param theta 0 < theta < pi - @param m - @return \int_0^\theta \sin^m x dx - **/ - if ((theta < 0) || (theta > M_PI)) - GSL_ERROR("Incomplete Wallis integral only allows theta in [0,pi]", GSL_EDOM); - if (theta < M_PI_2) - return 0.5 * beta_inc_unnormalized((m+1)/2.0, 0.5, gsl_pow_2(sin(theta))); - else - return 0.5 * (gsl_sf_beta((m+1)/2.0, 0.5) + beta_inc_unnormalized(0.5, (m+1)/2.0, gsl_pow_2(cos(theta)))); -} - -double planar_angle_to_solid_angle (double planar_angle, unsigned int dimension) -{ - return 2 * pow(M_PI, (dimension - 1)/2.0) - * incomplete_wallis_integral(planar_angle, dimension - 2) - / gsl_sf_gamma((dimension - 1)/2.0); -} - -double solid_angle_to_planar_angle (double solid_angle, unsigned int dimension) -{ - double f (double planar_angle, void* params) - { - return planar_angle_to_solid_angle(planar_angle, dimension) - solid_angle; - } - - gsl_function gsl_f = {&f}; - /* if (fabs(GSL_FN_EVAL(&gsl_f, 0))/surface_area_of_ball(dimension) < BISECTION_FUNCTION_EPSABS) */ - /* return 0; */ - /* else if (fabs(GSL_FN_EVAL(&gsl_f, M_PI))/surface_area_of_ball(dimension) < BISECTION_FUNCTION_EPSABS) */ - /* return M_PI; */ - /* else return bisection(&gsl_f, 0, M_PI); */ - if (solid_angle == 0) - return 0; - else if (solid_angle == surface_area_of_ball(dimension)) - return M_PI; - else return bisection(&gsl_f, 0, M_PI); -} diff --git a/src/nd-random.sc b/src/nd-random.sc new file mode 100644 index 0000000..3d96113 --- /dev/null +++ b/src/nd-random.sc @@ -0,0 +1,98 @@ +(pre-include "math.h") +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "gsl/gsl_sf_gamma.h") +(pre-include "gsl/gsl_vector.h") +(pre-include "utils.h") + +(define (random-direction-vector r x) + (void (const gsl-rng*) gsl-vector*) + "Generate a random vector distributed uniformly on the unit +sphere. Write the result to X." + (gsl-ran-dir-nd r (: x size) (: x data))) + +(define (rotate-from-nth-canonical x orient) + ((static void) gsl-vector* (const gsl-vector*)) + (let* ((n size-t (: x size)) + (xn double (gsl-vector-get x (- n 1))) + (mun double (gsl-vector-get orient (- n 1))) + (orient-sub gsl-vector-const-view + (gsl-vector-const-subvector orient 0 (- n 1))) + (b double (gsl-blas-dnrm2 (address-of (struct-get orient-sub vector)))) + (a double (/ (- (dot-product orient x) + (* xn mun)) + b)) + (s double (sqrt (- 1 (gsl-pow-2 mun))))) + (gsl-blas-daxpy (/ (+ (* xn s) + (* a (- mun 1))) + b) + orient + x) + (gsl-vector-set x + (- n 1) + (+ (gsl-vector-get x (- n 1)) + (* xn (- mun 1)) + (- (* a s)) + (- (/ (* mun (+ (* xn s) + (* a (- mun 1)))) + b)))))) + +(define (beta-inc-unnormalized a b x) ((static double) double double double) + (return (* (gsl-sf-beta-inc a b x) + (gsl-sf-beta a b)))) + +(define (incomplete-wallis-integral theta m) ((static double) double (unsigned int)) + "Return the incomplete Wallis integral \\int_0^\\theta \\sin^m x +dx. THETA should be in [0,pi]." + (cond + ((or (< theta 0) (> theta M-PI)) + (GSL-ERROR "Incomplete Wallis integral only allows theta in [0,pi]" GSL-EDOM)) + ((< theta M-PI-2) + (return (/ (beta-inc-unnormalized (* 0.5 (+ m 1)) 0.5 (gsl-pow-2 (sin theta))) + 2))) + (else + (return (/ (+ (gsl-sf-beta (* 0.5 (+ m 1)) 0.5) + (beta-inc-unnormalized 0.5 (* 0.5 (+ m 1)) (gsl-pow-2 (cos theta)))) + 2))))) + +(define (planar-angle->solid-angle planar-angle dimension) ((static double) double (unsigned int)) + (return (/ (* 2 + (pow M-PI (* 0.5 (- dimension 1))) + (incomplete-wallis-integral planar-angle (- dimension 2))) + (gsl-sf-gamma (* 0.5 (- dimension 1)))))) + +(define (solid-angle->planar-angle solid-angle dimension) ((static double) double (unsigned int)) + (define (f planar-angle params) (double double void*) + (return (- (planar-angle->solid-angle planar-angle dimension) + solid-angle))) + + (declare gsl-f (struct-variable gsl-function (function (address-of f)))) + ;; TODO: Equality comparisons to floating point values may be + ;; problematic. Fix it. + (cond + ((= solid-angle 0) (return 0)) + ((= solid-angle (surface-area-of-ball dimension)) (return M-PI)) + (else (return (bisection (address-of gsl-f) 0 M-PI))))) + +;; TODO: There is an edge case when mean is the (n-1)th canonical +;; basis vector. Fix it. +(define (hollow-cone-random-vector r mean theta-min theta-max x) + (void (const gsl-rng*) (const gsl-vector*) double double gsl-vector*) + ;; Generate random vector around the nth canonical basis vector. + (let* ((n size-t (: x size))) + (gsl-ran-dir-nd r (- n 1) (: x data)) + (gsl-vector-scale x (* (cos theta-max) + (tan (solid-angle->planar-angle + (gsl-ran-flat r + (planar-angle->solid-angle theta-min n) + (planar-angle->solid-angle theta-max n)) + n)))) + (gsl-vector-set x (- n 1) (cos theta-max))) + (gsl-vector-scale x (/ 1 (gsl-blas-dnrm2 x))) + + ;; Rotate to arbitrary basis. + (rotate-from-nth-canonical x mean)) + +(define (subsampling-random-vector r mean theta-max x) + (void (const gsl-rng*) (const gsl-vector*) double gsl-vector*) + (hollow-cone-random-vector r mean 0 theta-max x)) diff --git a/src/oracles.c b/src/oracles.c deleted file mode 100644 index 2b440eb..0000000 --- a/src/oracles.c +++ /dev/null @@ -1,108 +0,0 @@ -#include <math.h> -#include <gsl/gsl_math.h> -#include <gsl/gsl_randist.h> -#include "oracles.h" -#include "utils.h" - -double bernoulli_extent_generator -(const gsl_rng* r, double p, double r0, double r1) -{ - return gsl_ran_bernoulli(r, p) ? r1 : r0; -} - -double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension) -{ - return volume_of_ball(dimension) - * (p * gsl_pow_uint(r1, dimension) - + (1 - p) * gsl_pow_uint(r0, dimension)); -} - -double uniform_extent_generator (const gsl_rng* r, double a, double b) -{ - return gsl_ran_flat(r, a, b); -} - -// TODO: Verify the accuracy of this function for non-trivial a, b -double uniform_true_volume (double a, double b, unsigned int dimension) -{ - return exp(ln_volume_of_ball(dimension) + dimension*log(b) - log(dimension + 1)) - - exp(ln_volume_of_ball(dimension) + dimension*log(a) - log(dimension + 1)); -} - -double beta_extent_generator -(const gsl_rng* r, double alpha, double beta) -{ - return gsl_ran_beta(r, alpha, beta); -} - -double beta_true_volume (double alpha, double beta, unsigned int dimension) -{ - double vol = volume_of_ball(dimension); - for (unsigned int r=0; r<dimension; r++) - vol *= (alpha + r) / (alpha + beta + r); - return vol; -} - -static double infinity_norm (const gsl_vector* x) -{ - double max = fabs(gsl_vector_get(x, 0)); - for (int i=1; i<x->size; i++) - max = GSL_MAX(max, fabs(gsl_vector_get(x, i))); - return max; -} - -double cube_extent_oracle (const gsl_vector* x, double edge) -{ - return edge / 2 / infinity_norm(x); -} - -double cube_extent_oracle_with_center -(const gsl_vector* x, const gsl_vector* center, double edge) -{ - double min = (0.5*edge - GSL_SIGN(gsl_vector_get(x, 0))*gsl_vector_get(center, 0)) - / fabs(gsl_vector_get(x, 0)); - for (int i=1; i<center->size; i++) - min = GSL_MIN(min, (0.5*edge - GSL_SIGN(gsl_vector_get(x, i))*gsl_vector_get(center, i)) - / fabs(gsl_vector_get(x, i))); - return min; -} - -double cube_true_volume (double edge, unsigned int dimension) -{ - return gsl_pow_uint(edge, dimension); -} - -double cube_maximum_extent (double edge, unsigned int dimension) -{ - return edge * sqrt(dimension) / 2; -} - -double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes) -{ - double k = 0; - for (int i=0; i<axes->size; i++) - k += gsl_pow_2(gsl_vector_get(x, i) / gsl_vector_get(axes, i)); - return 1/sqrt(k); -} - -double ellipsoid_true_volume (const gsl_vector* axes) -{ - unsigned int dimension = axes->size; - double vol = volume_of_ball(dimension); - for (int i=0; i<dimension; i++) - vol *= gsl_vector_get(axes, i); - return vol; -} - -double spheroid_extent_oracle (const gsl_vector* x, double eccentricity) -{ - unsigned int dimension = x->size; - gsl_vector_const_view xsub = gsl_vector_const_subvector(x, 1, dimension - 1); - return 1/sqrt(gsl_pow_2(gsl_blas_dnrm2(&xsub.vector)) - + gsl_pow_2(gsl_vector_get(x, 0) / eccentricity)); -} - -double spheroid_true_volume (double eccentricity, unsigned int dimension) -{ - return volume_of_ball(dimension) * eccentricity; -} diff --git a/src/oracles.sc b/src/oracles.sc new file mode 100644 index 0000000..9680479 --- /dev/null +++ b/src/oracles.sc @@ -0,0 +1,90 @@ +(sc-include "macros/macros") + +(pre-include "math.h") +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "utils.h") + +(define (bernoulli-extent-generator r p r0 r1) (double (const gsl-rng*) double double double) + (return (if* (gsl-ran-bernoulli r p) r1 r0))) + +(define (bernoulli-true-volume p r0 r1 dimension) (double double double double (unsigned int)) + (return (* (volume-of-ball dimension) + (+ (* p (gsl-pow-uint r1 dimension)) + (* (- 1 p) (gsl-pow-uint r0 dimension)))))) + +(define (uniform-extent-generator r a b) (double (const gsl-rng*) double double) + (return (gsl-ran-flat r a b))) + +;; TODO: Verify the accuracy of this function for non-trivial a, b. +(define (uniform-true-volume a b dimension) (double double double (unsigned int)) + (return (- (exp (+ (ln-volume-of-ball dimension) + (* dimension (log b)) + (- (log (+ dimension 1))))) + (exp (+ (ln-volume-of-ball dimension) + (* dimension (log a)) + (- (log (+ dimension 1)))))))) + +(define (beta-extent-generator r alpha beta) (double (const gsl-rng*) double double) + (return (gsl-ran-beta r alpha beta))) + +(define (beta-true-volume alpha beta dimension) (double double double (unsigned int)) + (let* ((vol double (volume-of-ball dimension))) + (for-i r dimension + (set* vol (/ (+ alpha r) + (+ alpha beta r)))) + (return vol))) + +(define (infinity-norm x) ((static double) (const gsl-vector*)) + (let* ((max double (fabs (gsl-vector-get x 0)))) + ;; TODO: Start this loop from i = 1, not i = 0. That would be + ;; slightly faster. + (for-i i (: x size) + (set max (GSL-MAX max (fabs (gsl-vector-get x i))))) + (return max))) + +(define (cube-extent-oracle x edge) (double (const gsl-vector*) double) + (return (/ edge 2 (infinity-norm x)))) + +(sc-define-syntax (compute-cube-extent-oracle-minimizand i) + (/ (- (/ edge 2) + (* (GSL-SIGN (gsl-vector-get x i)) + (gsl-vector-get center i))) + (fabs (gsl-vector-get x i)))) + +(define (cube-extent-oracle-with-center x center edge) (double (const gsl-vector*) (const gsl-vector*) double) + (let* ((min double (compute-cube-extent-oracle-minimizand 0))) + ;; TODO: Start this loop from i = 1, not i = 0. That would be + ;; slightly faster. + (for-i i (: center size) + (set min (GSL-MIN min (compute-cube-extent-oracle-minimizand i)))) + (return min))) + +(define (cube-true-volume edge dimension) (double double (unsigned int)) + (return (gsl-pow-uint edge dimension))) + +(define (cube-maximum-extent edge dimension) (double double (unsigned int)) + (return (/ (* edge (sqrt dimension)) + 2))) + +(define (ellipsoid-extent-oracle x axes) (double (const gsl-vector*) (const gsl-vector*)) + (let* ((k double 0)) + (for-i i (: axes size) + (set+ k (gsl-pow-2 (/ (gsl-vector-get x i) + (gsl-vector-get axes i))))) + (return (/ (sqrt k))))) + +(define (ellipsoid-true-volume axes) (double (const gsl-vector*)) + (let* ((vol double (volume-of-ball (: axes size)))) + (for-i i (: axes size) + (set* vol (gsl-vector-get axes i))) + (return vol))) + +(define (spheroid-extent-oracle x eccentricity) (double (const gsl-vector*) double) + (let* ((xsub gsl-vector-const-view + (gsl-vector-const-subvector x 1 (- (: x size) 1)))) + (return (/ (sqrt (+ (gsl-pow-2 (gsl-blas-dnrm2 (address-of (struct-get xsub vector)))) + (gsl-pow-2 (/ (gsl-vector-get x 0) eccentricity)))))))) + +(define (spheroid-true-volume eccentricity dimension) (double double (unsigned int)) + (return (* (volume-of-ball dimension) eccentricity))) diff --git a/src/utils.c b/src/utils.c deleted file mode 100644 index b22e07b..0000000 --- a/src/utils.c +++ /dev/null @@ -1,160 +0,0 @@ -#include <gsl/gsl_blas.h> -#include <gsl/gsl_math.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_sf.h> -#include <gsl/gsl_sf_gamma.h> -#include "utils.h" - -// This function exists so that guile code can access the value of -// M_PI as provided by <gsl_math.h>. -double pi () -{ - /** - @return Value of pi - **/ - return M_PI; -} - -double volume_of_ball (unsigned int dimension) -{ - /** - @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.) - @return Volume of the unit ball - **/ - return exp(ln_volume_of_ball(dimension)); -} - -double ln_volume_of_ball (unsigned int dimension) -{ - /** - @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.) - @return Natural logarithm of the volume of the unit ball - **/ - return (dimension/2.0)*M_LNPI - gsl_sf_lngamma(1 + dimension/2.0); -} - -double surface_area_of_ball (unsigned int dimension) -{ - /** - @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.) - @return Surface area of the unit ball - **/ - return dimension * volume_of_ball(dimension); -} - -double ln_surface_area_of_ball (unsigned int dimension) -{ - /** - @param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.) - @return Natural logarithm of the surface area of the unit ball - **/ - return log(dimension) + ln_volume_of_ball(dimension); -} - -double lower_incomplete_gamma (double s, double x) -{ - return gsl_sf_gamma(s) * gsl_sf_gamma_inc_P(s, x); -} - -double angle_between_vectors (const gsl_vector* x, const gsl_vector* y) -{ - /** - @param x Vector 1 - @param y Vector 2 - @return Angle between the two vectors in the range [0, pi] - **/ - double dot_product; - gsl_blas_ddot(x, y, &dot_product); - // TODO: Is this a valid floating point comparison? - if (dot_product != 0) - return acos(dot_product / gsl_blas_dnrm2(x) / gsl_blas_dnrm2(y)); - else return 0; -} - -double dot_product (const gsl_vector* x, const gsl_vector* y) -{ - /** - @param x Vector 1 - @param y Vector 2 - @return Dot product between the two vectors - **/ - double result; - gsl_blas_ddot(x, y, &result); - return result; -} - -double gaussian_pdf (double x) -{ - /** - @param x - @return exp(-x^2/2) / sqrt(2*pi) - **/ - return gsl_ran_gaussian_pdf(x, 1.0); -} - -double gaussian_cdf (double x) -{ - /** - @param x - @return int_{-inf}^x gaussian_pdf(t) dt - **/ - return 0.5*(1 + gsl_sf_erf(x/M_SQRT2)); -} - -double rerror (double approx, double exact) -{ - /** - @param approx Approximate value - @param exact Exact value - @return Relative error - **/ - return fabs(1 - approx/exact); -} - -#define BISECTION_EPSABS 0 -#define BISECTION_EPSREL 1e-6 -#define BISECTION_MAX_ITERATIONS 1000 -#define BISECTION_FUNCTION_EPSABS 1e-9 -double bisection (gsl_function* f, double a, double b) -{ - const gsl_root_fsolver_type* T = gsl_root_fsolver_bisection; - gsl_root_fsolver* s = gsl_root_fsolver_alloc(T); - double solution; - int status; - - void error_handler (const char* reason, const char* file, int line, int gsl_errno) - { - fprintf(stderr, "Bisection error handler invoked.\n"); - fprintf(stderr, "f(%g) = %g\n", a, GSL_FN_EVAL(f, a)); - fprintf(stderr, "f(%g) = %g\n", b, GSL_FN_EVAL(f, b)); - fprintf(stderr, "gsl: %s:%d: ERROR: %s\n", file, line, reason); - abort(); - } - - gsl_error_handler_t* old_handler = gsl_set_error_handler(error_handler); - gsl_root_fsolver_set(s, f, a, b); - do { - status = gsl_root_fsolver_iterate (s); - solution = gsl_root_fsolver_root (s); - a = gsl_root_fsolver_x_lower (s); - b = gsl_root_fsolver_x_upper (s); - status = gsl_root_test_interval (a, b, BISECTION_EPSABS, BISECTION_EPSREL); - } while (status == GSL_CONTINUE); - gsl_root_fsolver_free(s); - gsl_set_error_handler(old_handler); - return solution; -} - -double bisection_rlimit (gsl_function* f, double a, double b) -{ - int sign = SIGNUM(GSL_FN_EVAL(f, a)); - for (int i=0; i<BISECTION_MAX_ITERATIONS; i++) - if (sign*GSL_FN_EVAL(f, b) > 0) b *= 2; - else return b; - fprintf(stderr, "Bisection bracketing failed\n"); - exit(1); -} -#undef BISECTION_EPSABS -#undef BISECTION_EPSREL -#undef BISECTION_MAX_ITERATIONS -#undef BISECTION_FUNCTION_EPSABS diff --git a/src/utils.sc b/src/utils.sc new file mode 100644 index 0000000..751cd77 --- /dev/null +++ b/src/utils.sc @@ -0,0 +1,123 @@ +(sc-include "macros/macros") + +(pre-include "gsl/gsl_math.h") +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "gsl/gsl_roots.h") +(pre-include "gsl/gsl_sf.h") + +(pre-define (SIGNUM x) + (if* (< x 0) -1 1)) + +;; This function exists so that guile code can access the value of +;; M_PI as provided by <gsl_math.h>. +(define (pi) (double) + "Return the value of pi." + (return M-PI)) + +(define (ln-volume-of-ball dimension) (double (unsigned int)) + "Return the natural logarithm of the volume of the unit +ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a +sphere, etc." + (return (- (* 0.5 dimension M-LNPI) + (gsl-sf-lngamma (+ 1 (* 0.5 dimension)))))) + +(define (volume-of-ball dimension) (double (unsigned int)) + "Return the volume of the unit ball of DIMENSION. DIMENSION of 2 +corresponds to a circle, 3 corresponds to a sphere, etc." + (return (exp (ln-volume-of-ball dimension)))) + +(define (ln-surface-area-of-ball dimension) (double (unsigned int)) + "Return the natural logarithm of the surface area of the unit +ball. DIMENSION of 2 corresponds to a circle, 3 corresponds to a +sphere, etc." + (return (+ (log dimension) + (ln-volume-of-ball dimension)))) + +(define (surface-area-of-ball dimension) (double (unsigned int)) + "Return the surface area of the unit ball of DIMENSION. DIMENSION of +2 corresponds to a circle, 3 corresponds to a sphere, etc." + (return (* dimension (volume-of-ball dimension)))) + +(define (lower-incomplete-gamma s x) (double double double) + (return (* (gsl-sf-gamma s) + (gsl-sf-gamma-inc-P s x)))) + +(define (angle-between-vectors x y) (double (const gsl-vector*) (const gsl-vector*)) + "Return the angle between vectors X and Y. The returned value is in +the range [0,pi]." + (declare dot-product double) + (gsl-blas-ddot x y (address-of dot-product)) + ;; TODO: Is this a valid floating point comparison? + (return (if* (= dot_product 0) + 0 + (acos (/ dot-product + (gsl-blas-dnrm2 x) + (gsl-blas-dnrm2 y)))))) + +(define (dot-product x y) (double (const gsl-vector*) (const gsl-vector*)) + "Return the dot product of vectors X and Y." + (declare result double) + (gsl-blas-ddot x y (address-of result)) + (return result)) + +(define (gaussian-pdf x) (double double) + "Return exp(-x^2/2) / sqrt(2*pi)" + (return (gsl-ran-gaussian-pdf x 1))) + +(define (gaussian-cdf x) (double double) + "Return \\int_{-\\inf}^x gaussian_pdf(t) dt." + (return (* 0.5 (+ 1 (gsl-sf-erf (/ x M-SQRT2)))))) + +(define (rerror approx exact) (double double double) + "Return the relative error between approximate value APPROX and +exact value EXACT." + (return (fabs (- 1 (/ approx exact))))) + +(sc-define-syntax* (with-root-fsolver solver solver-type function a b body ...) + (let ((solver-type (sc-gensym))) + `(begin + (let* ((,solver-type (const gsl-root-fsolver-type*) ,solver-type)) + (with-alloc solver gsl-root-fsolver* + (gsl-root-fsolver-alloc ,solver-type) + gsl-root-fsolver-free + (gsl-root-fsolver-set solver ,function ,a ,b) + ,@body))))) + +(sc-define-syntax* (with-error-handler handler body ...) + (let ((old-handler (sc-gensym))) + `(begin + (let* ((,old-handler gsl-error-handler-t* (gsl-set-error-handler ,handler))) + ,@body + (gsl-set-error-handler ,old-handler))))) + +(pre-let* (BISECTION-EPSABS 0 BISECTION-EPSREL 1e-6 BISECTION-MAX-ITERATIONS 1000) + (define (bisection f a b) (double gsl-function* double double) + (declare solution double) + (define (error-handler reason file line gsl-errno) (void (const char*) (const char*) int int) + (fprintf stderr "Bisection error handler invoked.\n") + (fprintf stderr "f(%g) = %g\n" a (GSL_FN_EVAL f a)) + (fprintf stderr "f(%g) = %g\n" b (GSL_FN_EVAL f b)) + (fprintf stderr "gsl: %s:%d: ERROR: %s\n" file line reason) + (abort)) + + (with-root-fsolver solver gsl-root-fsolver-bisection f a b + (with-error-handler error-handler + (do-while (= (gsl-root-test-interval (gsl-root-fsolver-x-lower solver) + (gsl-root-fsolver-x-upper solver) + BISECTION-EPSABS + BISECTION-EPSREL) + GSL-CONTINUE) + (gsl-root-fsolver-iterate solver))) + (set solution (gsl-root-fsolver-root solver))) + (return solution)) + + (define (bisection-rlimit f a b) (double gsl-function* double double) + (let* ((sign int (SIGNUM (GSL_FN_EVAL f a)))) + (for-i i BISECTION-MAX-ITERATIONS + (cond + ((> (* sign (GSL_FN_EVAL f b)) 0) + (set* b 2)) + (else (return b))))) + (fprintf stderr "Bisection bracketing failed\n") + (abort))) |