aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile15
-rw-r--r--extent-sampling.c161
-rw-r--r--extent-sampling.h39
-rw-r--r--nd-random.c101
-rw-r--r--nd-random.h15
-rw-r--r--oracles.c141
-rw-r--r--oracles.h37
-rw-r--r--utils.c181
-rw-r--r--utils.h28
9 files changed, 718 insertions, 0 deletions
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..81c9d98
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,15 @@
+CC = gcc
+CFLAGS = -Wall -D_GNU_SOURCE -g -fPIC -Werror
+LDFLAGS = -lm -lgsl -lgslcblas
+
+libextentsampling.so: extent-sampling.o nd-random.o gaussian-nd-random.o oracles.o utils.o
+ $(CC) $(CFLAGS) $(LDFLAGS) -shared -o $@ $^
+
+depends: $(wildcard *.c)
+ $(CC) $(CFLAGS) -MM $^ > depends
+
+include depends
+
+.PHONY: clean
+clean:
+ rm -vf *.o *.so depends
diff --git a/extent-sampling.c b/extent-sampling.c
new file mode 100644
index 0000000..eff3b3f
--- /dev/null
+++ b/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/extent-sampling.h b/extent-sampling.h
new file mode 100644
index 0000000..8ed3fb2
--- /dev/null
+++ b/extent-sampling.h
@@ -0,0 +1,39 @@
+#ifndef EXTENT_SAMPLING_H
+#define EXTENT_SAMPLING_H
+
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_rstat.h>
+#include <gsl/gsl_vector.h>
+
+typedef double (*extent_oracle_t) (const gsl_vector*);
+typedef double (*integrand_t) (double, const gsl_vector*);
+
+void init_random (void);
+
+double volume
+(extent_oracle_t extent_oracle, double true_volume,
+ const gsl_rng* r, unsigned int dimension, double rtol,
+ gsl_rstat_workspace* stats);
+
+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);
+
+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);
+
+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);
+
+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);
+
+#endif
diff --git a/nd-random.c b/nd-random.c
new file mode 100644
index 0000000..25abc7a
--- /dev/null
+++ b/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/nd-random.h b/nd-random.h
new file mode 100644
index 0000000..9637a0e
--- /dev/null
+++ b/nd-random.h
@@ -0,0 +1,15 @@
+#ifndef ND_RANDOM_H
+#define ND_RANDOM_H
+
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_vector.h>
+
+void random_direction_vector (const gsl_rng* r, gsl_vector* x);
+void hollow_cone_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_min, double theta_max, gsl_vector* x);
+void subsampling_random_vector (const gsl_rng* r, const gsl_vector* mean, double theta_max, gsl_vector* x);
+
+double planar_angle_to_solid_angle (double planar_angle, unsigned int dimension);
+double solid_angle_to_planar_angle (double solid_angle, unsigned int dimension);
+
+#endif
diff --git a/oracles.c b/oracles.c
new file mode 100644
index 0000000..80f401f
--- /dev/null
+++ b/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/oracles.h b/oracles.h
new file mode 100644
index 0000000..a661c94
--- /dev/null
+++ b/oracles.h
@@ -0,0 +1,37 @@
+#ifndef ORACLES_H
+#define ORACLES_H
+
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_rng.h>
+
+double bernoulli_extent_generator (const gsl_rng* r, double p, double r0, double r1);
+double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension);
+
+double uniform_extent_generator (const gsl_rng* r, double a, double b);
+double uniform_true_volume (double a, double b, unsigned int dimension);
+
+double beta_extent_generator (const gsl_rng* r, double alpha, double beta);
+double beta_true_volume (double alpha, double beta, unsigned int dimension);
+
+double symmetric_spiral_extent_oracle (const gsl_vector* x);
+
+double right_triangle_extent_oracle (const gsl_vector* x, double base, double height);
+double right_triangle_true_volume (double base, double height);
+
+double sphere_extent_oracle (const gsl_vector* x, double radius);
+double sphere_maximum_extent (double radius);
+
+double plane_extent_oracle (const gsl_vector* x, double displacement);
+
+double cube_extent_oracle (const gsl_vector* x, double edge);
+double cube_extent_oracle_with_center (const gsl_vector* x, const gsl_vector* center, double edge);
+double cube_true_volume (double edge, unsigned int dimension);
+double cube_maximum_extent (double edge, unsigned int dimension);
+
+double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes);
+double ellipsoid_true_volume (const gsl_vector* axes);
+
+double spheroid_extent_oracle (const gsl_vector* x, double eccentricity);
+double spheroid_true_volume (double eccentricity, unsigned int dimension);
+
+#endif
diff --git a/utils.c b/utils.c
new file mode 100644
index 0000000..0e0205f
--- /dev/null
+++ b/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
diff --git a/utils.h b/utils.h
new file mode 100644
index 0000000..220320e
--- /dev/null
+++ b/utils.h
@@ -0,0 +1,28 @@
+#ifndef UTILS_H
+#define UTILS_H
+
+#include <gsl/gsl_roots.h>
+#include <gsl/gsl_vector.h>
+
+#define SIGNUM(x) ((x) < 0 ? -1 : 1)
+
+double volume_of_ball (unsigned int dimension);
+double ln_volume_of_ball (unsigned int dimension);
+double surface_area_of_ball (unsigned int dimension);
+double ln_surface_area_of_ball (unsigned int dimension);
+double lower_incomplete_gamma (double s, double x);
+
+double angle_between_vectors (const gsl_vector* x, const gsl_vector* y);
+double dot_product (const gsl_vector* x, const gsl_vector* y);
+double gaussian_pdf (double x);
+double gaussian_cdf (double x);
+
+double rerror (double approx, double exact);
+
+double gprate (double start, double last, unsigned int n);
+double gpterm (double a, double r, int n);
+
+double bisection (gsl_function* f, double a, double b);
+double bisection_rlimit (gsl_function* f, double a, double b);
+
+#endif