blob: 71051fd3b61511e3cbaf3a11b095261cf4f7710b (
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
|
;;; 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 "stddef.h")
(pre-include "gsl/gsl_integration.h")
(pre-include "gsl/gsl_rstat.h")
(pre-include "gsl/gsl_vector.h")
(pre-include "nd-random.h")
(pre-include "utils.h")
(pre-include "extent-sampling.h")
(sc-define-syntax (with-integration-workspace var intervals body ...)
(with-alloc var gsl-integration-workspace*
(gsl-integration-workspace-alloc intervals)
gsl-integration-workspace-free
body ...))
(sc-define-syntax (invoke-extent-oracle extent-oracle r x)
((: extent-oracle oracle) r x (: extent-oracle params)))
(pre-define CONFIDENCE-INTERVAL-FACTOR 1.96)
(define (volume extent-oracle true-volume r dimension rtol stats)
(void (const extent-oracle-t*) double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
(let* ((vn double (ln-volume-of-ball dimension)))
(with-vector x dimension
(do-while (> (rerror (gsl-rstat-mean stats) true-volume)
rtol)
(random-direction-vector r x)
(gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
stats)))))
(sc-define-syntax (invoke-integrand integrand r x)
((: integrand integrand) r x (: integrand params)))
(pre-let* (INTEGRATION-INTERVALS 1000)
(define (integral-per-direction integrand direction n radius rtol)
(double (const integrand-t*) (const gsl-vector*) (unsigned int) double double)
(define (f r params) (double double void*)
(return (* (gsl-pow-uint r (- n 1))
(invoke-integrand integrand r direction))))
(declare gsl-f (struct-variable gsl-function (function (address-of f)))
result double
error double)
(with-integration-workspace w INTEGRATION-INTERVALS
(gsl-integration-qag (address-of gsl-f)
0
radius
0
rtol
INTEGRATION-INTERVALS
GSL-INTEG-GAUSS15
w
(address-of result)
(address-of error)))
(return (* (surface-area-of-ball n)
result))))
(define (integral integrand extent-oracle true-integral r dimension rtol stats)
(void (const integrand-t*) (const extent-oracle-t*) double (const gsl-rng*) (unsigned int) double gsl-rstat-workspace*)
(with-vector x dimension
(do-while (> (rerror (gsl-rstat-mean stats) true-integral)
rtol)
(random-direction-vector r x)
(gsl-rstat-add (integral-per-direction integrand x dimension
(invoke-extent-oracle extent-oracle r x) rtol)
stats))))
(define (volume-cone extent-oracle r mean omega-min omega-max number-of-samples variance)
(double (const extent-oracle-t*) (const gsl-rng*) (const gsl-vector*) double double (unsigned int) double*)
(declare volume double)
(let* ((dimension (unsigned int) (: mean size))
(vn double (ln-volume-of-ball dimension))
(theta-min double (solid-angle-fraction->planar-angle omega-min dimension))
(theta-max double (solid-angle-fraction->planar-angle omega-max dimension)))
(with-rstats stats
(with-vector x dimension
(for-i i number-of-samples
(hollow-cone-random-vector r mean theta-min theta-max x)
(gsl-rstat-add (exp (+ vn (* dimension (log (invoke-extent-oracle extent-oracle r x)))))
stats)))
(cond
(variance (set (pointer-get variance) (gsl-rstat-variance stats))))
(set volume (* (gsl-rstat-mean stats)
(- omega-max omega-min)))))
(return volume))
(define (volume-importance extent-oracle r mean samples-per-cone solid-angle-factor
solid-angle-threshold-exponent-factor number-of-samples)
(double (const extent-oracle-t*) (const gsl-rng*) (const gsl-vector*)
(unsigned int) double double (unsigned int*))
(define
volume double 0
omega-max double 0.5
N int 0)
(let* ((dimension int (: mean size)))
(do-while (> omega-max (pow 2 (- (* dimension solid-angle-threshold-exponent-factor))))
(let* ((omega-min double (/ omega-max solid-angle-factor)))
(set+ volume (volume-cone extent-oracle r mean omega-min omega-max samples-per-cone NULL))
(set+ N samples-per-cone)
(set omega-max omega-min))))
(set+ volume (volume-cone extent-oracle r mean 0 omega-max samples-per-cone NULL))
(set+ N samples-per-cone)
(cond
(number-of-samples (set (pointer-get number-of-samples) N)))
(return volume))
|
