diff options
Diffstat (limited to 'extent-sampling.c')
| -rw-r--r-- | extent-sampling.c | 161 |
1 files changed, 0 insertions, 161 deletions
diff --git a/extent-sampling.c b/extent-sampling.c deleted file mode 100644 index eff3b3f..0000000 --- a/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; -} |