diff options
author | Arun Isaac | 2021-02-05 15:00:27 +0530 |
---|---|---|
committer | Arun Isaac | 2021-02-05 15:00:27 +0530 |
commit | c31a57815be4d3ef5f6fd28cc9a34b3be55bfe56 (patch) | |
tree | c1be91a98018535ea5fd934e1986c3d3e5ba0128 /src/gaussian-nd-random.c | |
parent | 7352aa8150d0012d10be7b1aaa341be7c459a05a (diff) | |
download | nsmc-c31a57815be4d3ef5f6fd28cc9a34b3be55bfe56.tar.gz nsmc-c31a57815be4d3ef5f6fd28cc9a34b3be55bfe56.tar.lz nsmc-c31a57815be4d3ef5f6fd28cc9a34b3be55bfe56.zip |
Migrate C source to SC.
sph-sc is a scheme-like S-expression syntax for C. It elements much of
the pain and repetition involved in writing C syntax.
* src/extent-sampling.c, src/gaussian-nd-random.c, src/nd-random.c,
src/oracles.c, src/utils.c: Delete files.
* src/extent-sampling.sc, src/gaussian-nd-random.sc,
src/macros/macros.sc, src/nd-random.sc, src/oracles.sc, src/utils.sc:
New files.
* CMakeLists.txt: Generate C source files from SC source files.
Diffstat (limited to 'src/gaussian-nd-random.c')
-rw-r--r-- | src/gaussian-nd-random.c | 51 |
1 files changed, 0 insertions, 51 deletions
diff --git a/src/gaussian-nd-random.c b/src/gaussian-nd-random.c deleted file mode 100644 index 16abc94..0000000 --- a/src/gaussian-nd-random.c +++ /dev/null @@ -1,51 +0,0 @@ -#include <gsl/gsl_blas.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_sf_gamma.h> -#include "nd-random.h" -#include "utils.h" - -#define INTEGRATION_EPSABS 1e-3 -#define INTEGRATION_EPSREL 1e-3 - -unsigned int shifted_gaussian_random_vector -(const gsl_rng* r, const gsl_vector* mean, double theta_max, double standard_deviation, gsl_vector* x) -{ - for (int i=0; i<x->size; i++) - gsl_vector_set(x, i, gsl_ran_gaussian(r, standard_deviation)); - gsl_vector_add(x, mean); - gsl_blas_dscal(1/gsl_blas_dnrm2(x), x); - if (angle_between_vectors(x, mean) > theta_max) - return 1 + shifted_gaussian_random_vector(r, mean, theta_max, standard_deviation, x); - else return 1; -} - -double shifted_gaussian_pdf (double theta, double mean, double theta_max, - double standard_deviation, unsigned int dimension, gsl_integration_workspace* w) -{ - // NOTE: Due to the rejection sampling, the returned result is only - // proportional to the actual PDF. - if (theta > theta_max) return 0; - - double projection = mean*cos(theta); - double integrand (double r, void* params) - { - return gsl_pow_uint(r, dimension - 1) - * gaussian_pdf((r - projection) / standard_deviation); - } - - gsl_function gsl_integrand = {&integrand}; - double result, error; - /* gsl_integration_qagiu(&gsl_integrand, 0, INTEGRATION_EPSABS, INTEGRATION_EPSREL, */ - /* w->limit, w, &result, &error); */ - double rmax = 0.5 * (projection + sqrt(gsl_pow_2(projection) + 4*(dimension-1)*gsl_pow_2(standard_deviation))); - double halfwidth = 3 * standard_deviation; - gsl_integration_qag(&gsl_integrand, - rmax - halfwidth < 0 ? 0 : rmax - halfwidth, - rmax + halfwidth, - INTEGRATION_EPSABS, INTEGRATION_EPSREL, - w->limit, GSL_INTEG_GAUSS41, w, &result, &error); - return result - * gaussian_pdf(sqrt(gsl_pow_2(mean) - gsl_pow_2(projection)) / standard_deviation) - / gsl_pow_uint(M_SQRT2*sqrt(M_PI), dimension - 2) - / gsl_pow_uint(standard_deviation, dimension); -} |