aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--experiments/cube.c46
-rw-r--r--experiments/integral.c113
-rw-r--r--experiments/spheroid.c57
-rw-r--r--experiments/volume.c106
4 files changed, 322 insertions, 0 deletions
diff --git a/experiments/cube.c b/experiments/cube.c
new file mode 100644
index 0000000..845688c
--- /dev/null
+++ b/experiments/cube.c
@@ -0,0 +1,46 @@
+/* How does the cube fare with unbiased extent sampling?
+ */
+
+#include <stdio.h>
+#include <gsl/gsl_vector.h>
+#include "oracles.h"
+#include "extent-sampling.h"
+
+#define DIMENSION_START 5
+#define DIMENSION_STEP 5
+#define DIMENSION_STOP 50
+
+#define TRIALS 100
+
+#define RTOL 0.1
+#define EDGE 1.0
+
+double cube_oracle (const gsl_vector* x) {
+ return cube_extent_oracle(x, EDGE);
+}
+
+int main ()
+{
+ gsl_rng_env_setup();
+ gsl_rng* r = gsl_rng_alloc(gsl_rng_default);
+
+ FILE* fp = fopen("cube-window.dat", "w");
+ printf("dimension\tsamples\n");
+ fprintf(fp, "dimension\tsamples\n");
+ for (int dim=DIMENSION_START; dim<=DIMENSION_STOP; dim+=DIMENSION_STEP) {
+ unsigned int total_samples=0;
+ for (int i=0; i<TRIALS; i++) {
+ unsigned int number_of_samples;
+ volume_window(cube_oracle, cube_true_volume(EDGE, dim), r, dim, RTOL, &number_of_samples);
+ total_samples += number_of_samples;
+ }
+ double average_samples = ((double)total_samples)/TRIALS;
+ printf("%d\t%g\n", dim, average_samples);
+ fprintf(fp, "%d\t%g\n", dim, average_samples);
+ fflush(NULL);
+ }
+ fclose(fp);
+
+ gsl_rng_free(r);
+ return 0;
+}
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;
+}
diff --git a/experiments/spheroid.c b/experiments/spheroid.c
new file mode 100644
index 0000000..c2c8403
--- /dev/null
+++ b/experiments/spheroid.c
@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <math.h>
+#include <time.h>
+#include "extent-sampling.h"
+#include "oracles.h"
+
+#define DIMENSION_START 5
+#define DIMENSION_STEP 5
+#define DIMENSION_STOP 100
+
+#define RTOL 0.1
+#define SAMPLES_PER_CONE 1000
+#define SOLID_ANGLE_FACTOR 1.1
+
+void run_experiment
+(double eccentricity, double solid_angle_threshold_exponent_factor, FILE* fp)
+{
+ gsl_rng_env_setup();
+ gsl_rng* r = gsl_rng_alloc(gsl_rng_default);
+
+ spheroid_params params = {eccentricity};
+ extent_oracle_t oracle = {spheroid_extent_oracle, &params};
+
+ fprintf(fp, "# spheroid eccentricity = %g\n", eccentricity);
+ fprintf(fp, "# samples per cone = %d\n", SAMPLES_PER_CONE);
+ fprintf(fp, "# solid angle threshold = 2^(-%gn)\n", solid_angle_threshold_exponent_factor);
+ fprintf(fp, "# solid angle factor = %g\n", SOLID_ANGLE_FACTOR);
+ fprintf(fp, "dimension\tvolume\tsamples\ttime\n");
+ for (int dim=DIMENSION_START; dim<=DIMENSION_STOP; dim+=DIMENSION_STEP) {
+ gsl_vector* mean = gsl_vector_alloc(dim);
+ gsl_vector_set_basis(mean, 0);
+
+ unsigned int samples;
+ clock_t begin = clock();
+ double volume = volume_importance(&oracle, r, mean, SAMPLES_PER_CONE,
+ SOLID_ANGLE_FACTOR, solid_angle_threshold_exponent_factor, &samples);
+ clock_t end = clock();
+ fprintf(fp, "%d\t%g\t%d\t%g\n", dim, 2*volume / spheroid_true_volume(dim, &params),
+ samples, (double)(end-begin)/CLOCKS_PER_SEC);
+ fflush(NULL);
+
+ gsl_vector_free(mean);
+ }
+ gsl_rng_free(r);
+}
+
+int main ()
+{
+ FILE* fp = fopen("spheroid-eccentricity10.dat", "w");
+ run_experiment(10, 2, fp);
+ fclose(fp);
+
+ fp = fopen("spheroid-eccentricity100.dat", "w");
+ run_experiment(100, 6, fp);
+ fclose(fp);
+ return 0;
+}
diff --git a/experiments/volume.c b/experiments/volume.c
new file mode 100644
index 0000000..3a20d43
--- /dev/null
+++ b/experiments/volume.c
@@ -0,0 +1,106 @@
+#include <math.h>
+#include <stdio.h>
+
+#include <gsl/gsl_randist.h>
+
+#include "extent-sampling.h"
+#include "oracles.h"
+#include "utils.h"
+
+#define RTOLS 3
+#define RTOL1 0.2
+#define RTOL2 0.1
+#define RTOL3 0.05
+
+typedef double (*true_volume_t) (unsigned int);
+gsl_rng* r;
+
+#define MAX_EXTENT 1.0
+double uniform_oracle (const gsl_vector* x) {
+ return gsl_ran_flat(r, 0, MAX_EXTENT);
+}
+
+double uniform_true_volume_dim (unsigned int dim)
+{
+ return uniform_true_volume(0, MAX_EXTENT, dim);
+}
+#undef MAX_EXTENT
+
+#define ALPHA 2
+#define BETA 2
+double beta_oracle (const gsl_vector* x) {
+ return beta_extent_generator(r, ALPHA, BETA);
+}
+
+double beta_true_volume_dim (unsigned int dim)
+{
+ return beta_true_volume(ALPHA, BETA, dim);
+}
+#undef ALPHA
+#undef BETA
+
+double arcsine_oracle (const gsl_vector* x) {
+ return beta_extent_generator(r, 0.5, 0.5);
+}
+
+double arcsine_true_volume_dim (unsigned int dim)
+{
+ return beta_true_volume(0.5, 0.5, dim);
+}
+
+#define EDGE 1.0
+double cube_oracle (const gsl_vector* x) {
+ return cube_extent_oracle(x, EDGE);
+}
+
+double cube_true_volume_dim (unsigned int dim)
+{
+ return cube_true_volume(EDGE, dim);
+}
+#undef EDGE
+
+void run_experiment (extent_oracle_t oracle, true_volume_t true_volume, const char* filename_prefix,
+ int dim_start, int dim_step, int dim_stop)
+{
+ 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=dim_start; dim<=dim_stop; dim+=dim_step) {
+ printf("%d\t", dim);
+ for (int i=0; i<RTOLS; i++) {
+ volume(oracle, true_volume(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("# uniform\n");
+ run_experiment(uniform_oracle, uniform_true_volume_dim, "volume-uniform", 5, 5, 100);
+ printf("# beta\n");
+ run_experiment(beta_oracle, beta_true_volume_dim, "volume-beta", 5, 5, 100);
+ printf("# arcsine\n");
+ run_experiment(arcsine_oracle, arcsine_true_volume_dim, "volume-arcsine", 5, 5, 100);
+ printf("# cube\n");
+ run_experiment(cube_oracle, cube_true_volume_dim, "volume-cube", 10, 10, 50);
+
+ gsl_rng_free(r);
+ return 0;
+}
generated by cgit v1.2.3 (git 2.45.2) at 2024-12-22 09:29:23 +0000