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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
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;
}
|
