blob: 689ac9d4e8d20ba95882320e54ffae1c288f424e (
plain)
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
134
135
136
137
138
139
140
141
142
143
144
145
146
|
;;; nsmc --- n-sphere Monte Carlo method
;;; Copyright © 2021 Arun I <arunisaac@systemreboot.net>
;;; Copyright © 2021 Murugesan Venkatapathi <murugesh@iisc.ac.in>
;;;
;;; This file is part of nsmc.
;;;
;;; nsmc is free software: you can redistribute it and/or modify it
;;; under the terms of the GNU General Public License as published by
;;; the Free Software Foundation, either version 3 of the License, or
;;; (at your option) any later version.
;;;
;;; nsmc is distributed in the hope that it will be useful, but
;;; WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
;;; General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with nsmc. If not, see <https://www.gnu.org/licenses/>.
(sc-include "macros/macros")
(pre-include "math.h")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_poly.h")
(pre-include "gsl/gsl_randist.h")
(pre-include "oracles.h")
(pre-include "utils.h")
(define (bernoulli-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const bernoulli-params*) (convert-type -params bernoulli-params*)))
(return (if* (gsl-ran-bernoulli r (: params p))
(: params r1)
(: params r0)))))
(define (bernoulli-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params (const bernoulli-params*) (convert-type -params bernoulli-params*)))
(return (* (volume-of-ball dimension)
(+ (* (: params p) (gsl-pow-uint (: params r1) dimension))
(* (- 1 (: params p)) (gsl-pow-uint (: params r0) dimension)))))))
(define (uniform-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const uniform-params*) (convert-type -params uniform-params*)))
(return (gsl-ran-flat r (: params a) (: params b)))))
;; TODO: Verify the accuracy of this function for non-trivial a, b.
(define (uniform-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params uniform-params* (convert-type -params uniform-params*)))
(return (* (volume-of-ball dimension)
(- (/ (gsl-pow-uint (: params b) (+ dimension 1))
(+ dimension 1))
(/ (gsl-pow-uint (: params a) (+ dimension 1))
(+ dimension 1)))))))
(define (beta-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const beta-params*) (convert-type -params beta-params*)))
(return (gsl-ran-beta r (: params alpha) (: params beta)))))
(define (beta-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params (const beta-params*) (convert-type -params beta-params*))
(vol double (volume-of-ball dimension)))
(for-i r dimension
(set* vol (/ (+ (: params alpha) r)
(+ (: params alpha) (: params beta) r))))
(return vol)))
(define (infinity-norm x) ((static double) (const gsl-vector*))
(let* ((max double (fabs (gsl-vector-get x 0))))
;; TODO: Start this loop from i = 1, not i = 0. That would be
;; slightly faster.
(for-i i (: x size)
(set max (GSL-MAX max (fabs (gsl-vector-get x i)))))
(return max)))
(define (cube-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const cube-params*) (convert-type -params cube-params*)))
(return (/ (: params edge) 2 (infinity-norm x)))))
(sc-define-syntax (compute-cube-extent-oracle-minimizand i)
(/ (- (/ (: params edge) 2)
(* (GSL-SIGN (gsl-vector-get x i))
(gsl-vector-get (: params center) i)))
(fabs (gsl-vector-get x i))))
(define (cube-extent-oracle-with-center r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const cube-params*) (convert-type -params cube-params*))
(min double (compute-cube-extent-oracle-minimizand 0)))
;; TODO: Start this loop from i = 1, not i = 0. That would be
;; slightly faster.
(for-i i (: (: params center) size)
(set min (GSL-MIN min (compute-cube-extent-oracle-minimizand i))))
(return min)))
(define (cube-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params (const cube-params*) (convert-type -params cube-params*)))
(return (gsl-pow-uint (: params edge) dimension))))
(define (ellipsoid-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*))
(k double 0))
(for-i i (: (: params axes) size)
(set+ k (gsl-pow-2 (/ (gsl-vector-get x i)
(gsl-vector-get (: params axes) i)))))
(return (/ (sqrt k)))))
(define (ellipsoid-extent-oracle-with-center r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*))
;; a, b and c are the coefficients of a quadratic equation
;; ax^2 + bx + c = 0.
(a double 0)
(b double 0)
(c double 0))
(for-i i (: (: params axes) size)
(set+ a (gsl-pow-2 (/ (gsl-vector-get x i)
(gsl-vector-get (: params axes) i))))
(set+ b (/ (* (gsl-vector-get x i)
(gsl-vector-get (: params center) i))
(gsl-pow-2 (gsl-vector-get (: params axes) i))))
(set+ c (gsl-pow-2 (/ (gsl-vector-get (: params center) i)
(gsl-vector-get (: params axes) i)))))
(set* b 2)
(set- c 1)
(declare k1 double
k2 double)
(gsl-poly-solve-quadratic a b c (address-of k1) (address-of k2))
;; k2 > k1, k1 < 0, k2 > 0. k2 is our desired root.
(return k2)))
(define (ellipsoid-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params (const ellipsoid-params*) (convert-type -params ellipsoid-params*))
(vol double (volume-of-ball (: (: params axes) size))))
(for-i i (: (: params axes) size)
(set* vol (gsl-vector-get (: params axes) i)))
(return vol)))
(define (spheroid-extent-oracle r x -params) (double (const gsl-rng*) (const gsl-vector*) void*)
(let* ((params (const spheroid-params*) (convert-type -params spheroid-params*))
(xsub gsl-vector-const-view
(gsl-vector-const-subvector x 1 (- (: x size) 1))))
(return (/ (sqrt (+ (gsl-pow-2 (gsl-blas-dnrm2 (address-of (struct-get xsub vector))))
(gsl-pow-2 (/ (gsl-vector-get x 0)
(: params eccentricity)))))))))
(define (spheroid-true-volume dimension -params) (double (unsigned int) void*)
(let* ((params (const spheroid-params*) (convert-type -params spheroid-params*)))
(return (* (volume-of-ball dimension)
(: params eccentricity)))))
|
