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.

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;
+}