diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/extent-sampling.c | 161 | ||||
| -rw-r--r-- | src/gaussian-nd-random.c | 51 | ||||
| -rw-r--r-- | src/nd-random.c | 101 | ||||
| -rw-r--r-- | src/oracles.c | 141 | ||||
| -rw-r--r-- | src/utils.c | 181 |
5 files changed, 635 insertions, 0 deletions
diff --git a/src/extent-sampling.c b/src/extent-sampling.c new file mode 100644 index 0000000..eff3b3f --- /dev/null +++ b/src/extent-sampling.c @@ -0,0 +1,161 @@ +/* + 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/gaussian-nd-random.c b/src/gaussian-nd-random.c new file mode 100644 index 0000000..16abc94 --- /dev/null +++ b/src/gaussian-nd-random.c @@ -0,0 +1,51 @@ +#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/nd-random.c b/src/nd-random.c new file mode 100644 index 0000000..25abc7a --- /dev/null +++ b/src/nd-random.c @@ -0,0 +1,101 @@ +#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/oracles.c b/src/oracles.c new file mode 100644 index 0000000..80f401f --- /dev/null +++ b/src/oracles.c @@ -0,0 +1,141 @@ +#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; +} + +double symmetric_spiral_extent_oracle (const gsl_vector* x) +{ + double angle = atan2(gsl_vector_get(x, 1), gsl_vector_get(x, 0)); + return fabs(angle); +} + +double right_triangle_extent_oracle (const gsl_vector* x, double base, double height) +{ + double max_angle = atan(height / base); + double angle = atan2(gsl_vector_get(x, 1), gsl_vector_get(x, 0)); + return (angle > 0) && (angle < max_angle) ? base / cos(angle) : 0; +} + +double right_triangle_true_volume (double base, double height) +{ + return 0.5 * base * height; +} + +double sphere_extent_oracle (const gsl_vector* x, double radius) +{ + return radius; +} + +double sphere_maximum_extent (double radius) +{ + return radius; +} + +double plane_extent_oracle (const gsl_vector* x, double displacement) +{ + return displacement / fabs(gsl_vector_get(x, 0)); +} + +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/utils.c b/src/utils.c new file mode 100644 index 0000000..0e0205f --- /dev/null +++ b/src/utils.c @@ -0,0 +1,181 @@ +#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); +} + +double gprate (double start, double last, unsigned int n) +{ + /** + @param start First (0th) term of the geometric progression + @param last Last ('N-1'th) term of the geometric progression + @param n Number of terms in the geometric progression + @return The multiplicative factor of the geometric progression + **/ + return pow(last/start, 1.0/(n - 1)); +} + +double gpterm (double a, double r, int n) +{ + /** + @param a First (0th) term of the geometric progression + @param r Factor of the geometric progression + @return The Nth term of the geometric progression, with N starting from 0 + **/ + return a * gsl_pow_uint(r, n); +} + +#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 |