aboutsummaryrefslogtreecommitdiff
path: root/utils.c
diff options
context:
space:
mode:
authorArun Isaac2021-02-03 19:40:10 +0530
committerArun Isaac2021-02-03 19:40:10 +0530
commit983f091af7849f259366f641b438ac1fee3df940 (patch)
treefd42196b6031638fd906afd8acbc5c26e34ce305 /utils.c
parent707afdfe66dd9c931af77e4ff73186fd89ed27fd (diff)
downloadnsmc-983f091af7849f259366f641b438ac1fee3df940.tar.gz
nsmc-983f091af7849f259366f641b438ac1fee3df940.tar.lz
nsmc-983f091af7849f259366f641b438ac1fee3df940.zip
Move source files and headers to separate directories.
* extent-sampling.h, gaussian-nd-random.h, nd-random.h, oracles.h, utils.h: Move into include directory. * extent-sampling.c, gaussian-nd-random.c, nd-random.c, oracles.c, utils.c: Move into src directory. * CMakeLists.txt: Set include as include directory. Look for source files inside src directory.
Diffstat (limited to 'utils.c')
-rw-r--r--utils.c181
1 files changed, 0 insertions, 181 deletions
diff --git a/utils.c b/utils.c
deleted file mode 100644
index 0e0205f..0000000
--- a/utils.c
+++ /dev/null
@@ -1,181 +0,0 @@
-#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