1 files changed, 113 insertions, 0 deletions
diff --git a/experiments/integral.c b/experiments/integral.c
new file mode 100644
index 0000000..46827d9
--- /dev/null
+++ b/experiments/integral.c
@@ -0,0 +1,113 @@
+#include <complex.h>
+#include <math.h>
+#include <omp.h>
+#include <string.h>
+
+#include <gsl/gsl_sf.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_poly.h>
+#include <gsl/gsl_randist.h>
+
+#include "extent-sampling.h"
+#include "utils.h"
+
+#define DIMENSION_START 5
+#define DIMENSION_STEP 5
+#define DIMENSION_STOP 100
+
+#define RTOLS 3
+#define RTOL1 0.2
+#define RTOL2 0.1
+#define RTOL3 0.05
+
+#define MAX_EXTENT 1.0
+
+#define POLYNOMIAL_DEGREE 3
+double poly_coeffs[POLYNOMIAL_DEGREE + 1] = {-0.09375, 0.6875, -1.5, 1};
+
+typedef double (*true_integral_t) (unsigned int);
+gsl_rng* r;
+
+double uniform_oracle (const gsl_vector* x) {
+  return gsl_ran_flat(r, 0, MAX_EXTENT);
+}
+
+double polynomial_integrand  (double r, const gsl_vector* x) {
+  return gsl_poly_eval(poly_coeffs, POLYNOMIAL_DEGREE + 1, r);
+}
+
+double polynomial_true_integral (unsigned int n) {
+  double integrand[POLYNOMIAL_DEGREE + 1];
+  for (int k=0; k<=POLYNOMIAL_DEGREE; k++)
+    integrand[k] = poly_coeffs[k]/(n + k)/(n + k + 1);
+  return surface_area_of_ball(n)
+    / MAX_EXTENT
+    * gsl_pow_uint(MAX_EXTENT, n + 1)
+    * gsl_poly_eval(integrand, POLYNOMIAL_DEGREE+1, MAX_EXTENT);
+}
+
+double gaussian_integrand  (double r, const gsl_vector* x) {
+  return exp(-gsl_pow_2(r)/2);
+}
+
+double gaussian_true_integral (unsigned int n) {
+  return pow(2, 0.5*n - 1)
+    * surface_area_of_ball(n)
+    / MAX_EXTENT
+    * (MAX_EXTENT*lower_incomplete_gamma(0.5*n, 0.5*gsl_pow_2(MAX_EXTENT))
+       - sqrt(2) * lower_incomplete_gamma(0.5*(n + 1), 0.5*gsl_pow_2(MAX_EXTENT)));
+}
+
+double x_coordinate_integrand (double r, const gsl_vector* x) {
+  return fabs(r*gsl_vector_get(x, 0));
+}
+
+double x_coordinate_true_integral (unsigned int n) {
+  return surface_area_of_ball(n+3) / 2 / (n+2) / gsl_pow_2(M_PI)
+    * gsl_pow_uint(MAX_EXTENT, n+1)
+    / MAX_EXTENT;
+}
+
+void run_experiment (integrand_t integrand, extent_oracle_t oracle, true_integral_t true_integral, const char* filename_prefix)
+{
+  gsl_rstat_workspace* stats = gsl_rstat_alloc();
+  printf("dimension\trtol=%g\trtol=%g\trtol=%g\n", RTOL1, RTOL2, RTOL3);
+  double rtol[RTOLS] = {RTOL1, RTOL2, RTOL3};
+  FILE *fp[RTOLS];
+  for (int i=0; i<RTOLS; i++) {
+    char *filename;
+    asprintf(&filename, "%s-%g.dat", filename_prefix, rtol[i]);
+    fp[i] = fopen(filename, "w");
+    fprintf(fp[i], "dimension\tsamples\n");
+    free(filename);
+  }
+  for (int dim=DIMENSION_START; dim<=DIMENSION_STOP; dim+=DIMENSION_STEP) {
+    printf("%d\t", dim);
+    for (int i=0; i<RTOLS; i++) {
+      integral(integrand, oracle, true_integral(dim), r, dim, rtol[i], stats);
+      fprintf(fp[i], "%d\t%zd\n", dim, gsl_rstat_n(stats));
+      printf("%zd\t", gsl_rstat_n(stats));
+    }
+    printf("\n");
+    gsl_rstat_reset(stats);
+  }
+  for (int i=0; i<RTOLS; i++) fclose(fp[i]);
+  gsl_rstat_free(stats);
+}
+
+int main (int argc, char** argv)
+{
+  gsl_rng_env_setup();
+  r = gsl_rng_alloc(gsl_rng_default);
+
+  printf("# polynomial\n");
+  run_experiment(polynomial_integrand, uniform_oracle, polynomial_true_integral, "integral-polynomial");
+  printf("# gaussian\n");
+  run_experiment(gaussian_integrand, uniform_oracle, gaussian_true_integral, "integral-gaussian");
+  printf("# x-coordinate\n");
+  run_experiment(x_coordinate_integrand, uniform_oracle, x_coordinate_true_integral, "integral-x-coordinate");
+
+  gsl_rng_free(r);
+  return 0;
+}