diff options
Diffstat (limited to 'src/utils.c')
| -rw-r--r-- | src/utils.c | 181 |
1 files changed, 181 insertions, 0 deletions
diff --git a/src/utils.c b/src/utils.c new file mode 100644 index 0000000..0e0205f --- /dev/null +++ b/src/utils.c @@ -0,0 +1,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 |