1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
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;
}
|