1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
;;; SPDX-FileCopyrightText: 2020 Arvydas Silanskas
;;; SPDX-FileCopyrightText: 2020 Linas Vepštas
;;; SPDX-FileCopyrightText: 2024 Bradley J Lucier
;;; SPDX-License-Identifier: MIT
;;;
;;; sphere.scm
;;; Uniform distributions on a sphere, and a ball.
;;; Submitted for inclusion in srfi-194
;;;
;;; Algorithm based on BoxMeuller as described in
;;; http://extremelearning.com.au/how-to-generate-uniformly-random-points-on-n-spheres-and-n-balls/
;;;
;;; make-sphere-generator N - return a generator of points uniformly
;;; distributed on an N-dimensional sphere.
;;; This implements the BoxMeuller algorithm, that is, of normalizing
;;; N+1 Gaussian random variables.
(define (make-sphere-generator arg)
(cond
((and (integer? arg)
(exact? arg)
(positive? arg))
(make-ellipsoid-generator* (make-vector (+ 1 arg) 1.0)))
(else
(error "make-sphere-generator: The argument must be a positive exact integer: " arg))))
(define (make-ellipsoid-generator arg)
(define (return-error)
(error "make-ellipsoid-generator: The argument must be a vector of real numbers that are finite and positive when converted to inexact: " arg))
(if (and (vector? arg)
(vector-every real? arg))
(let ((inexact-arg (vector-map inexact arg)))
(if (vector-every (lambda (x)
(and (positive? x) (finite? x)))
inexact-arg)
(make-ellipsoid-generator* inexact-arg)
(return-error)))
(return-error)))
;;; -----------------------------------------------
;;; Generator of points uniformly distributed on an N-dimensional ellipsoid.
;;;
;;; The `axes` should be a vector of floats, specifying the axes of the
;;; ellipsoid. The algorithm used is an accept/reject sampling algo,
;;; wherein the acceptance rate is proportional to the measure of a
;;; surface element on the ellipsoid. The measure is straight-forward to
;;; arrive at, and the 3D case is described by `mercio` in detail at
;;; https://math.stackexchange.com/questions/973101/how-to-generate-points-uniformly-distributed-on-the-surface-of-an-ellipsoid
;;;
;;; Note that sampling means that complexity goes as
;;; O(B/A x C/A x D/A x ...) where `A` is the shorest axis,
;;; and `B`, `C`, `D`, ... are the other axes. Maximum performance
;;; is achieved on spheres, which is the case used in make-ball-generator
(define (make-ellipsoid-generator* axes)
(let ((gauss (make-normal-generator))
(uniform (make-random-real-generator 0. 1.)) ;; should really be from a separate stream
(min-axis (vector-fold min +inf.0 axes)))
(define (sphere)
(let* ((point
(vector-map (lambda (_) (gauss)) axes))
(norm-inverse
(/ (sqrt (vector-fold (lambda (sum x)
(+ sum (square x)))
0.
point)))))
(vector-map (lambda (x) (* x norm-inverse)) point)))
(define (ellipsoid-distance ray)
(sqrt (vector-fold
(lambda (sum x a) (+ sum (square (/ x a))))
0. ray axes)))
(define (keep point)
(< (uniform)
(* min-axis (ellipsoid-distance point))))
(define (sample)
(let ((point (sphere)))
(if (keep point)
point
(sample))))
(lambda ()
(vector-map * (sample) axes))))
;;;-----------------------------------------------
;;; make-ball-generator N - return a generator of points
;;; inside an N-dimensional ball. It's based on the algorithm in
;;;
;;; An Efficient Method for Generating Points
;;; Uniformly Distributed in Hyperellipsoids
;;; Jean Dezert and Christian Musso
;;; https://www.onera.fr/sites/default/files/297/C013_-_Dezert_-_YBSTributeMonterey2001.pdf
;;;
;;; which in turn is based on the Harman-Lacko-Voelker Dropped Coordinate method for
;;; generating points uniformly inside the unit ball in N dimensions.
(define (make-ball-generator arg)
(define (return-error)
(error "make-ball-generator: The argument must be either an exact, positive, integer or a vector of real numbers that are positive and finite when converted to inexact: " arg))
(if (and (integer? arg)
(exact? arg)
(positive? arg))
(make-ball-generator* (make-vector arg 1.0))
(if (and (vector? arg)
(vector-every real? arg))
(let ((inexact-arg (vector-map inexact arg)))
(if (vector-every (lambda (x)
(and (positive? x)
(finite? x)))
inexact-arg)
(make-ball-generator* inexact-arg)
(return-error)))
(return-error))))
(define (make-ball-generator* axes)
(let* ((sphere-generator
;; returns vectors with (vector-length axes) + 2 elements
(make-sphere-generator (+ (vector-length axes) 1))))
(lambda ()
;; returns vectors with (vector-length axes) elements
(vector-map (lambda (el axis)
(* el axis))
(sphere-generator)
axes))))
|
