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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_gamma.h>
#include "utils.h"
// This function exists so that guile code can access the value of
// M_PI as provided by <gsl_math.h>.
double pi ()
{
/**
@return Value of pi
**/
return M_PI;
}
double volume_of_ball (unsigned int dimension)
{
/**
@param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
@return Volume of the unit ball
**/
return exp(ln_volume_of_ball(dimension));
}
double ln_volume_of_ball (unsigned int dimension)
{
/**
@param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
@return Natural logarithm of the volume of the unit ball
**/
return (dimension/2.0)*M_LNPI - gsl_sf_lngamma(1 + dimension/2.0);
}
double surface_area_of_ball (unsigned int dimension)
{
/**
@param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
@return Surface area of the unit ball
**/
return dimension * volume_of_ball(dimension);
}
double ln_surface_area_of_ball (unsigned int dimension)
{
/**
@param dimension Dimension of the ball (2 for a circle, 3 for a sphere, etc.)
@return Natural logarithm of the surface area of the unit ball
**/
return log(dimension) + ln_volume_of_ball(dimension);
}
double lower_incomplete_gamma (double s, double x)
{
return gsl_sf_gamma(s) * gsl_sf_gamma_inc_P(s, x);
}
double angle_between_vectors (const gsl_vector* x, const gsl_vector* y)
{
/**
@param x Vector 1
@param y Vector 2
@return Angle between the two vectors in the range [0, pi]
**/
double dot_product;
gsl_blas_ddot(x, y, &dot_product);
// TODO: Is this a valid floating point comparison?
if (dot_product != 0)
return acos(dot_product / gsl_blas_dnrm2(x) / gsl_blas_dnrm2(y));
else return 0;
}
double dot_product (const gsl_vector* x, const gsl_vector* y)
{
/**
@param x Vector 1
@param y Vector 2
@return Dot product between the two vectors
**/
double result;
gsl_blas_ddot(x, y, &result);
return result;
}
double gaussian_pdf (double x)
{
/**
@param x
@return exp(-x^2/2) / sqrt(2*pi)
**/
return gsl_ran_gaussian_pdf(x, 1.0);
}
double gaussian_cdf (double x)
{
/**
@param x
@return int_{-inf}^x gaussian_pdf(t) dt
**/
return 0.5*(1 + gsl_sf_erf(x/M_SQRT2));
}
double rerror (double approx, double exact)
{
/**
@param approx Approximate value
@param exact Exact value
@return Relative error
**/
return fabs(1 - approx/exact);
}
double gprate (double start, double last, unsigned int n)
{
/**
@param start First (0th) term of the geometric progression
@param last Last ('N-1'th) term of the geometric progression
@param n Number of terms in the geometric progression
@return The multiplicative factor of the geometric progression
**/
return pow(last/start, 1.0/(n - 1));
}
double gpterm (double a, double r, int n)
{
/**
@param a First (0th) term of the geometric progression
@param r Factor of the geometric progression
@return The Nth term of the geometric progression, with N starting from 0
**/
return a * gsl_pow_uint(r, n);
}
#define BISECTION_EPSABS 0
#define BISECTION_EPSREL 1e-6
#define BISECTION_MAX_ITERATIONS 1000
#define BISECTION_FUNCTION_EPSABS 1e-9
double bisection (gsl_function* f, double a, double b)
{
const gsl_root_fsolver_type* T = gsl_root_fsolver_bisection;
gsl_root_fsolver* s = gsl_root_fsolver_alloc(T);
double solution;
int status;
void error_handler (const char* reason, const char* file, int line, int gsl_errno)
{
fprintf(stderr, "Bisection error handler invoked.\n");
fprintf(stderr, "f(%g) = %g\n", a, GSL_FN_EVAL(f, a));
fprintf(stderr, "f(%g) = %g\n", b, GSL_FN_EVAL(f, b));
fprintf(stderr, "gsl: %s:%d: ERROR: %s\n", file, line, reason);
abort();
}
gsl_error_handler_t* old_handler = gsl_set_error_handler(error_handler);
gsl_root_fsolver_set(s, f, a, b);
do {
status = gsl_root_fsolver_iterate (s);
solution = gsl_root_fsolver_root (s);
a = gsl_root_fsolver_x_lower (s);
b = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (a, b, BISECTION_EPSABS, BISECTION_EPSREL);
} while (status == GSL_CONTINUE);
gsl_root_fsolver_free(s);
gsl_set_error_handler(old_handler);
return solution;
}
double bisection_rlimit (gsl_function* f, double a, double b)
{
int sign = SIGNUM(GSL_FN_EVAL(f, a));
for (int i=0; i<BISECTION_MAX_ITERATIONS; i++)
if (sign*GSL_FN_EVAL(f, b) > 0) b *= 2;
else return b;
fprintf(stderr, "Bisection bracketing failed\n");
exit(1);
}
#undef BISECTION_EPSABS
#undef BISECTION_EPSREL
#undef BISECTION_MAX_ITERATIONS
#undef BISECTION_FUNCTION_EPSABS
|
