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.sc | |
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.sc')
-rw-r--r-- | src/gaussian-nd-random.sc | 72 |
1 files changed, 72 insertions, 0 deletions
diff --git a/src/gaussian-nd-random.sc b/src/gaussian-nd-random.sc new file mode 100644 index 0000000..bfde566 --- /dev/null +++ b/src/gaussian-nd-random.sc @@ -0,0 +1,72 @@ +(sc-include "macros/macros") + +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_integration.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "utils.h") + +(pre-define INTEGRATION-EPSABS 1e-3) +(pre-define INTEGRATION-EPSREL 1e-3) + +(define (shifted-gaussian-random-vector r mean theta-max standard-deviation x) + ((unsigned int) (const gsl-rng*) (const gsl-vector*) double double gsl-vector*) + "Return a random vector generated by the shifted Gaussian method." + (for-i i (: x size) + (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) + (cond + ((> (angle-between-vectors x mean) theta-max) + (return (+ 1 (shifted-gaussian-random-vector r mean theta-max standard-deviation x)))) + (else (return 1)))) + +(define (shifted-gaussian-pdf theta mean theta-max standard-deviation dimension w) + (double double double double double (unsigned int) gsl-integration-workspace*) + "Return the unnormalized probability density function at THETA for a +random vector generated using the shifted Gaussian method." + (cond + ((> theta theta-max) (return 0)) + (else + (define projection double (* mean (cos theta))) + (define (integrand r params) (double double void*) + (return (* (gsl-pow-uint r (- dimension 1)) + (gaussian-pdf (/ (- r projection) + standard-deviation))))) + + (declare gsl-integrand (struct-variable gsl-function (function (address-of integrand))) + result double + error double) + + ;; (gsl-integration-qagiu (address-of gsl-integrand) + ;; 0 + ;; INTEGRATION-EPSABS + ;; INTEGRATION-EPSREL + ;; (: w limit) + ;; w + ;; (address-of result) + ;; (address-of error)) + (let* ((rmax double (/ (+ projection + (sqrt (+ (gsl-pow-2 projection) + (* 4 + (- dimension 1) + (gsl-pow-2 standard-deviation))))) + 2)) + (half-width double (* 3 standard-deviation))) + (gsl-integration-qag (address-of gsl-integrand) + (if* (< (- rmax half-width) 0) + 0 + (- rmax half-width)) + (+ rmax half-width) + INTEGRATION-EPSABS + INTEGRATION-EPSREL + (: w limit) + GSL-INTEG-GAUSS41 + w + (address-of result) + (address-of error))) + (return (/ (* result + (gaussian-pdf (/ (sqrt (- (gsl-pow-2 mean) + (gsl-pow-2 projection))) + standard-deviation))) + (gsl-pow-uint (sqrt (* 2 M-PI)) (- dimension 2)) + (gsl-pow-uint standard-deviation dimension)))))) |