diff options
Diffstat (limited to 'src/extent-sampling.c')
| -rw-r--r-- | src/extent-sampling.c | 161 |
1 files changed, 161 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; +} |