diff options
Diffstat (limited to 'src/nd-random.sc')
| -rw-r--r-- | src/nd-random.sc | 98 |
1 files changed, 98 insertions, 0 deletions
diff --git a/src/nd-random.sc b/src/nd-random.sc new file mode 100644 index 0000000..3d96113 --- /dev/null +++ b/src/nd-random.sc @@ -0,0 +1,98 @@ +(pre-include "math.h") +(pre-include "gsl/gsl_blas.h") +(pre-include "gsl/gsl_randist.h") +(pre-include "gsl/gsl_sf_gamma.h") +(pre-include "gsl/gsl_vector.h") +(pre-include "utils.h") + +(define (random-direction-vector r x) + (void (const gsl-rng*) gsl-vector*) + "Generate a random vector distributed uniformly on the unit +sphere. Write the result to X." + (gsl-ran-dir-nd r (: x size) (: x data))) + +(define (rotate-from-nth-canonical x orient) + ((static void) gsl-vector* (const gsl-vector*)) + (let* ((n size-t (: x size)) + (xn double (gsl-vector-get x (- n 1))) + (mun double (gsl-vector-get orient (- n 1))) + (orient-sub gsl-vector-const-view + (gsl-vector-const-subvector orient 0 (- n 1))) + (b double (gsl-blas-dnrm2 (address-of (struct-get orient-sub vector)))) + (a double (/ (- (dot-product orient x) + (* xn mun)) + b)) + (s double (sqrt (- 1 (gsl-pow-2 mun))))) + (gsl-blas-daxpy (/ (+ (* xn s) + (* a (- mun 1))) + b) + orient + x) + (gsl-vector-set x + (- n 1) + (+ (gsl-vector-get x (- n 1)) + (* xn (- mun 1)) + (- (* a s)) + (- (/ (* mun (+ (* xn s) + (* a (- mun 1)))) + b)))))) + +(define (beta-inc-unnormalized a b x) ((static double) double double double) + (return (* (gsl-sf-beta-inc a b x) + (gsl-sf-beta a b)))) + +(define (incomplete-wallis-integral theta m) ((static double) double (unsigned int)) + "Return the incomplete Wallis integral \\int_0^\\theta \\sin^m x +dx. THETA should be in [0,pi]." + (cond + ((or (< theta 0) (> theta M-PI)) + (GSL-ERROR "Incomplete Wallis integral only allows theta in [0,pi]" GSL-EDOM)) + ((< theta M-PI-2) + (return (/ (beta-inc-unnormalized (* 0.5 (+ m 1)) 0.5 (gsl-pow-2 (sin theta))) + 2))) + (else + (return (/ (+ (gsl-sf-beta (* 0.5 (+ m 1)) 0.5) + (beta-inc-unnormalized 0.5 (* 0.5 (+ m 1)) (gsl-pow-2 (cos theta)))) + 2))))) + +(define (planar-angle->solid-angle planar-angle dimension) ((static double) double (unsigned int)) + (return (/ (* 2 + (pow M-PI (* 0.5 (- dimension 1))) + (incomplete-wallis-integral planar-angle (- dimension 2))) + (gsl-sf-gamma (* 0.5 (- dimension 1)))))) + +(define (solid-angle->planar-angle solid-angle dimension) ((static double) double (unsigned int)) + (define (f planar-angle params) (double double void*) + (return (- (planar-angle->solid-angle planar-angle dimension) + solid-angle))) + + (declare gsl-f (struct-variable gsl-function (function (address-of f)))) + ;; TODO: Equality comparisons to floating point values may be + ;; problematic. Fix it. + (cond + ((= solid-angle 0) (return 0)) + ((= solid-angle (surface-area-of-ball dimension)) (return M-PI)) + (else (return (bisection (address-of gsl-f) 0 M-PI))))) + +;; TODO: There is an edge case when mean is the (n-1)th canonical +;; basis vector. Fix it. +(define (hollow-cone-random-vector r mean theta-min theta-max x) + (void (const gsl-rng*) (const gsl-vector*) double double gsl-vector*) + ;; Generate random vector around the nth canonical basis vector. + (let* ((n size-t (: x size))) + (gsl-ran-dir-nd r (- n 1) (: x data)) + (gsl-vector-scale x (* (cos theta-max) + (tan (solid-angle->planar-angle + (gsl-ran-flat r + (planar-angle->solid-angle theta-min n) + (planar-angle->solid-angle theta-max n)) + n)))) + (gsl-vector-set x (- n 1) (cos theta-max))) + (gsl-vector-scale x (/ 1 (gsl-blas-dnrm2 x))) + + ;; Rotate to arbitrary basis. + (rotate-from-nth-canonical x mean)) + +(define (subsampling-random-vector r mean theta-max x) + (void (const gsl-rng*) (const gsl-vector*) double gsl-vector*) + (hollow-cone-random-vector r mean 0 theta-max x)) |