aboutsummaryrefslogtreecommitdiff
path: root/src/extent-sampling.c
diff options
context:
space:
mode:
authorArun Isaac2021-02-03 19:40:10 +0530
committerArun Isaac2021-02-03 19:40:10 +0530
commit983f091af7849f259366f641b438ac1fee3df940 (patch)
treefd42196b6031638fd906afd8acbc5c26e34ce305 /src/extent-sampling.c
parent707afdfe66dd9c931af77e4ff73186fd89ed27fd (diff)
downloadnsmc-983f091af7849f259366f641b438ac1fee3df940.tar.gz
nsmc-983f091af7849f259366f641b438ac1fee3df940.tar.lz
nsmc-983f091af7849f259366f641b438ac1fee3df940.zip
Move source files and headers to separate directories.
* extent-sampling.h, gaussian-nd-random.h, nd-random.h, oracles.h, utils.h: Move into include directory. * extent-sampling.c, gaussian-nd-random.c, nd-random.c, oracles.c, utils.c: Move into src directory. * CMakeLists.txt: Set include as include directory. Look for source files inside src directory.
Diffstat (limited to 'src/extent-sampling.c')
-rw-r--r--src/extent-sampling.c161
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;
+}