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
|
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include "oracles.h"
#include "utils.h"
double bernoulli_extent_generator
(const gsl_rng* r, double p, double r0, double r1)
{
return gsl_ran_bernoulli(r, p) ? r1 : r0;
}
double bernoulli_true_volume (double p, double r0, double r1, unsigned int dimension)
{
return volume_of_ball(dimension)
* (p * gsl_pow_uint(r1, dimension)
+ (1 - p) * gsl_pow_uint(r0, dimension));
}
double uniform_extent_generator (const gsl_rng* r, double a, double b)
{
return gsl_ran_flat(r, a, b);
}
// TODO: Verify the accuracy of this function for non-trivial a, b
double uniform_true_volume (double a, double b, unsigned int dimension)
{
return exp(ln_volume_of_ball(dimension) + dimension*log(b) - log(dimension + 1))
- exp(ln_volume_of_ball(dimension) + dimension*log(a) - log(dimension + 1));
}
double beta_extent_generator
(const gsl_rng* r, double alpha, double beta)
{
return gsl_ran_beta(r, alpha, beta);
}
double beta_true_volume (double alpha, double beta, unsigned int dimension)
{
double vol = volume_of_ball(dimension);
for (unsigned int r=0; r<dimension; r++)
vol *= (alpha + r) / (alpha + beta + r);
return vol;
}
double symmetric_spiral_extent_oracle (const gsl_vector* x)
{
double angle = atan2(gsl_vector_get(x, 1), gsl_vector_get(x, 0));
return fabs(angle);
}
double right_triangle_extent_oracle (const gsl_vector* x, double base, double height)
{
double max_angle = atan(height / base);
double angle = atan2(gsl_vector_get(x, 1), gsl_vector_get(x, 0));
return (angle > 0) && (angle < max_angle) ? base / cos(angle) : 0;
}
double right_triangle_true_volume (double base, double height)
{
return 0.5 * base * height;
}
double sphere_extent_oracle (const gsl_vector* x, double radius)
{
return radius;
}
double sphere_maximum_extent (double radius)
{
return radius;
}
double plane_extent_oracle (const gsl_vector* x, double displacement)
{
return displacement / fabs(gsl_vector_get(x, 0));
}
double infinity_norm (const gsl_vector* x)
{
double max = fabs(gsl_vector_get(x, 0));
for (int i=1; i<x->size; i++)
max = GSL_MAX(max, fabs(gsl_vector_get(x, i)));
return max;
}
double cube_extent_oracle (const gsl_vector* x, double edge)
{
return edge / 2 / infinity_norm(x);
}
double cube_extent_oracle_with_center
(const gsl_vector* x, const gsl_vector* center, double edge)
{
double min = (0.5*edge - GSL_SIGN(gsl_vector_get(x, 0))*gsl_vector_get(center, 0))
/ fabs(gsl_vector_get(x, 0));
for (int i=1; i<center->size; i++)
min = GSL_MIN(min, (0.5*edge - GSL_SIGN(gsl_vector_get(x, i))*gsl_vector_get(center, i))
/ fabs(gsl_vector_get(x, i)));
return min;
}
double cube_true_volume (double edge, unsigned int dimension)
{
return gsl_pow_uint(edge, dimension);
}
double cube_maximum_extent (double edge, unsigned int dimension)
{
return edge * sqrt(dimension) / 2;
}
double ellipsoid_extent_oracle (const gsl_vector* x, const gsl_vector* axes)
{
double k = 0;
for (int i=0; i<axes->size; i++)
k += gsl_pow_2(gsl_vector_get(x, i) / gsl_vector_get(axes, i));
return 1/sqrt(k);
}
double ellipsoid_true_volume (const gsl_vector* axes)
{
unsigned int dimension = axes->size;
double vol = volume_of_ball(dimension);
for (int i=0; i<dimension; i++)
vol *= gsl_vector_get(axes, i);
return vol;
}
double spheroid_extent_oracle (const gsl_vector* x, double eccentricity)
{
unsigned int dimension = x->size;
gsl_vector_const_view xsub = gsl_vector_const_subvector(x, 1, dimension - 1);
return 1/sqrt(gsl_pow_2(gsl_blas_dnrm2(&xsub.vector))
+ gsl_pow_2(gsl_vector_get(x, 0) / eccentricity));
}
double spheroid_true_volume (double eccentricity, unsigned int dimension)
{
return volume_of_ball(dimension) * eccentricity;
}
|
