Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementation of procedures in [pm_distGamma](@ref pm_distGamma).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the URNG.
28 : #if RNGD_ENABLED
29 : #define RNG
30 : #elif RNGF_ENABLED || RNGX_ENABLED
31 : #define RNG rng,
32 : #elif setGammaRand_ENABLED
33 : #error "Unrecognized interface."
34 : #endif
35 : !%%%%%%%%%%%%%%%%%%%%%%%
36 : #if getGammaLogPDFNF_ENABLED
37 : !%%%%%%%%%%%%%%%%%%%%%%%
38 :
39 17975924 : CHECK_ASSERTION(__LINE__, kappa > 0._RKC, SK_"@getGammaLogPDFNF(): The condition `kappa > 0.` must hold. kappa = "//getStr(kappa))
40 : #if KD_ENABLED
41 8988082 : logPDFNF = -log_gamma(kappa)
42 : #elif KS_ENABLED
43 8987842 : CHECK_ASSERTION(__LINE__, invSigma > 0._RKC, SK_"@getGammaLogPDFNF(): The condition invSigma > 0.` must hold. invSigma = "//getStr(invSigma))
44 8987842 : logPDFNF = getGammaLogPDFNF(kappa)
45 8987842 : if (invSigma /= 1._RKC) logPDFNF = logPDFNF + log(invSigma)
46 : #else
47 : #error "Unrecognized interface."
48 : #endif
49 :
50 : !%%%%%%%%%%%%%%%%%%%%%
51 : #elif getGammaLogPDF_ENABLED
52 : !%%%%%%%%%%%%%%%%%%%%%
53 :
54 : real(RKC) :: kappa_def, invSigma_def
55 2991515 : kappa_def = 1._RKC; if (present(kappa)) kappa_def = kappa
56 2991515 : invSigma_def = 1._RKC; if (present(invSigma)) invSigma_def = invSigma
57 2991515 : call setGammaLogPDF(logPDF, x, getGammaLogPDFNF(kappa_def, invSigma_def), kappa_def, invSigma_def)
58 :
59 : !%%%%%%%%%%%%%%%%%%%%%
60 : #elif setGammaLogPDF_ENABLED
61 : !%%%%%%%%%%%%%%%%%%%%%
62 :
63 : #if DDD_ENABLED
64 81 : CHECK_ASSERTION(__LINE__, x > 0._RKC, SK_"@setGammaLogPDF(): The condition `x > 0.` must hold. x = "//getStr(x))
65 81 : logPDF = -x
66 : #elif NKD_ENABLED
67 80 : CHECK_ASSERTION(__LINE__, x > 0._RKC, SK_"@setGammaLogPDF(): The condition `x > 0.` must hold. x = "//getStr(x))
68 80 : CHECK_ASSERTION(__LINE__, kappa > 0._RKC, SK_"@setGammaLogPDF(): The condition `kappa > 0.` must hold. kappa = "//getStr(kappa))
69 240 : CHECK_ASSERTION(__LINE__, abs(logPDFNF - getGammaLogPDFNF(kappa)) < sqrt(epsilon(0._RKC)), \
70 : SK_"@setGammaLogPDF(): The condition `abs(logPDFNF - getGammaLogPDFNF(kappa)) < sqrt(epsilon(0._RKC)` must hold. logPDFNF, getGammaLogPDFNF(kappa) = "// \
71 : getStr([logPDFNF, getGammaLogPDFNF(kappa)]))
72 80 : logPDF = logPDFNF + (kappa - 1._RKC) * log(x) - x
73 : #elif NKS_ENABLED
74 : real(RKC) :: y
75 2995611 : y = x * invSigma
76 2995611 : CHECK_ASSERTION(__LINE__, x > 0._RKC, SK_"@setGammaLogPDF(): The condition `x > 0.` must hold. x = "//getStr(x))
77 2995611 : CHECK_ASSERTION(__LINE__, kappa > 0._RKC, SK_"@setGammaLogPDF(): The condition `kappa > 0.` must hold. kappa = "//getStr(kappa))
78 2995611 : CHECK_ASSERTION(__LINE__, invSigma > 0._RKC, SK_"@setGammaLogPDF(): The condition `invSigma > 0.` must hold. invSigma = "//getStr(invSigma))
79 8986833 : CHECK_ASSERTION(__LINE__, abs(logPDFNF - getGammaLogPDFNF(kappa, invSigma)) < sqrt(epsilon(0._RKC)), \
80 : SK_"@setGammaLogPDF(): The condition `abs(logPDFNF - getGammaLogPDFNF(kappa, invSigma)) < sqrt(epsilon(0._RKC)` must hold. logPDFNF, getGammaLogPDFNF(kappa, invSigma) = "// \
81 : getStr([logPDFNF, getGammaLogPDFNF(kappa, invSigma)]))
82 2995611 : logPDF = logPDFNF + (kappa - 1._RKC) * log(y) - y
83 : #else
84 : #error "Unrecognized interface."
85 : #endif
86 :
87 : !%%%%%%%%%%%%%%%%%%
88 : #elif getGammaCDF_ENABLED
89 : !%%%%%%%%%%%%%%%%%%
90 :
91 : integer(IK) :: info
92 : real(RKC) :: xnormed
93 4017 : if (present(invSigma)) then
94 4017 : xnormed = x * invSigma
95 : else
96 0 : xnormed = x
97 : end if
98 4017 : if (present(kappa)) then
99 4017 : call setGammaCDF(cdf, xnormed, log_gamma(kappa), kappa, info)
100 : else
101 0 : call setGammaCDF(cdf, xnormed, info)
102 : end if
103 4017 : if (info < 0_IK) error stop MODULE_NAME//SK_"@getGammaCDF(): The computation of the regularized Lower Incomplete Gamma function failed. This can happen if `kappa` is too large."
104 :
105 : !%%%%%%%%%%%%%%%%%%
106 : #elif setGammaCDF_ENABLED
107 : !%%%%%%%%%%%%%%%%%%
108 :
109 : #if DD_ENABLED
110 : real(RKC), parameter :: kappa = 1._RKC, logGammaKappa = log_gamma(kappa)
111 1 : CHECK_ASSERTION(__LINE__, x > 0._RKC, SK_"@setGammaCDF(): The condition `x > 0.` must hold. x = "//getStr(x)) ! fpp
112 1 : call setGammaIncLow(cdf, x, logGammaKappa, kappa, info)
113 : #else
114 54229 : CHECK_ASSERTION(__LINE__, x > 0._RKC, SK_"@setGammaCDF(): The condition `x > 0.` must hold. x = "//getStr(x)) ! fpp
115 54229 : CHECK_ASSERTION(__LINE__, kappa > 0._RKC, SK_"@setGammaCDF(): The condition `kappa > 0.` must hold. kappa = "//getStr(kappa)) ! fpp
116 162687 : CHECK_ASSERTION(__LINE__, abs(log_gamma(kappa) - logGammaKappa) < 100 * epsilon(0._RKC), \
117 : SK_"@setGammaCDF(): The condition `abs(log_gamma(kappa) - logGammaKappa) < 100 * epsilon(0._RKC)` must hold. log_gamma(kappa), logGammaKappa = "//\
118 : getStr([log_gamma(kappa), logGammaKappa])) ! fpp
119 : #if KD_ENABLED
120 4017 : call setGammaIncLow(cdf, x, logGammaKappa, kappa, info)
121 : #elif KS_ENABLED
122 50212 : CHECK_ASSERTION(__LINE__, invSigma > 0._RKC, SK_"@setGammaCDF(): The condition `invSigma > 0.` must hold. invSigma = "//getStr(invSigma)) ! fpp
123 50212 : call setGammaIncLow(cdf, x * invSigma, logGammaKappa, kappa, info)
124 : #else
125 : #error "Unrecognized interface."
126 : #endif
127 : #endif
128 :
129 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 : #elif setGammaRand_ENABLED && (RNGD_ENABLED || RNGF_ENABLED || RNGX_ENABLED) && KR_ENABLED
131 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 :
133 : real(RKC) , parameter :: ONE_THIRD = 1._RKC / 3._RKC
134 : real(RKC) :: invSqrt9KappaMinusOneThird
135 : real(RKC) :: kappaMinusOneThird
136 : real(RKC) :: normrnd, unifrnd
137 1362 : CHECK_ASSERTION(__LINE__, 0._RKC < kappa .and. 0._RKC < sigma, SK_"@setGammaRand(): The condition `0 < kappa .and. 0 < sigma` must hold. kappa = "//getStr([kappa, sigma])) ! Must appear here.
138 : #define GET_PARAM(offset) \
139 : kappaMinusOneThird = kappa - ONE_THIRD + offset; invSqrt9KappaMinusOneThird = 1._RKC / sqrt(9 * kappaMinusOneThird);
140 454 : if (1._RKC <= kappa) then
141 440 : GET_PARAM(0._RKC) ! fpp
142 : #if D1_ENABLED
143 : #define GET_RAND(i) rand(i)
144 : block
145 : integer(IK) :: irand
146 40011 : do irand = 1_IK, size(rand, 1, IK)
147 : #elif D0_ENABLED
148 : #define GET_RAND(i) rand
149 : #else
150 : #error "Unrecognized interface."
151 : #endif
152 9 : do
153 : !call setRandKappaGe1()
154 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 : ! If only macros existed in Fortran...
156 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157 41131 : call setNormRand(RNG normrnd)
158 41131 : GET_RAND(irand) = 1._RKC + invSqrt9KappaMinusOneThird * normrnd
159 41131 : if (GET_RAND(irand) <= 0._RKC) cycle
160 41089 : call setUnifRand(RNG unifrnd); unifrnd = 1._RKC - unifrnd
161 41089 : if (unifrnd < 1._RKC - 0.0331_RKC * normrnd**4) then
162 37667 : GET_RAND(irand) = kappaMinusOneThird * sigma * GET_RAND(irand)**3
163 37667 : exit
164 : end if
165 3422 : GET_RAND(irand) = GET_RAND(irand)**3
166 3422 : if (log(unifrnd) < 0.5_RKC * normrnd**2 + kappaMinusOneThird * (1._RKC - GET_RAND(irand) + log(GET_RAND(irand)))) then
167 2775 : GET_RAND(irand) = kappaMinusOneThird * sigma * GET_RAND(irand)
168 2775 : exit
169 : end if
170 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171 : end do
172 : #if D1_ENABLED
173 : end do
174 : end block
175 : #endif
176 : else
177 2 : GET_PARAM(1._RKC) ! fpp
178 : #if D1_ENABLED
179 : block
180 : integer(IK) :: irand
181 5009 : do irand = 1_IK, size(rand, 1, IK)
182 : #endif
183 3 : do
184 : !call setRandKappaGe1()
185 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186 : ! If only macros existed in Fortran...
187 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
188 5120 : call setNormRand(RNG normrnd)
189 5120 : GET_RAND(irand) = 1._RKC + invSqrt9KappaMinusOneThird * normrnd
190 5120 : if (GET_RAND(irand) <= 0._RKC) cycle
191 5119 : call setUnifRand(RNG unifrnd); unifrnd = 1._RKC - unifrnd
192 5119 : if (unifrnd < 1._RKC - 0.0331_RKC * normrnd**4) then
193 4693 : GET_RAND(irand) = kappaMinusOneThird * sigma * GET_RAND(irand)**3
194 4693 : exit
195 : end if
196 426 : GET_RAND(irand) = GET_RAND(irand)**3
197 426 : if (log(unifrnd) < 0.5_RKC * normrnd**2 + kappaMinusOneThird * (1._RKC - GET_RAND(irand) + log(GET_RAND(irand)))) then
198 315 : GET_RAND(irand) = kappaMinusOneThird * sigma * GET_RAND(irand)
199 315 : exit
200 : end if
201 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202 111 : call setUnifRand(RNG unifrnd)
203 111 : unifrnd = (1._RKC - unifrnd)**(1._RKC / kappa)
204 112 : GET_RAND(irand) = GET_RAND(irand) * unifrnd
205 : end do
206 : #if D1_ENABLED
207 : end do
208 : end block
209 : #endif
210 : end if
211 :
212 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213 : #elif setGammaRand_ENABLED && (RNGD_ENABLED || RNGF_ENABLED || RNGX_ENABLED) && KI_ENABLED
214 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215 :
216 : ! Alternative algorithm when `kappa` is a positive integer.
217 : ! The corresponding interfaces are yet to be implemented.
218 : ! However, the benefit of this integer `kappa` is unclear.
219 : real(RKC) :: unifrnd(7)
220 : CHECK_ASSERTION(__LINE__, 0._RKC < kappa .and. 0._RKC < sigma, SK_"@setGammaRand(): The condition `0 < kappa .and. 0 < sigma` must hold. kappa = "//getStr([kappa, sigma]))
221 : if (kappa < 6_IK) then
222 : call setUnifRand(RNG unifrnd(1:kappa))
223 : rand = -log(product(unifrnd(1:kappa)))
224 : else ! use rejection sampling
225 : do
226 : call setUnifRand(RNG unifrnd(1:2))
227 : unifrnd(1:2) = 2 * unifrnd(1:2) - 1._RKC
228 : if (dot_product(unifrnd(1:2), unifrnd(1:2)) > 1._RKC) cycle
229 : unifrnd(3) = unifrnd(2) / unifrnd(1)
230 : unifrnd(4) = kappa - 1._RKC
231 : unifrnd(5) = sqrt(2 * unifrnd(4) + 1._RKC)
232 : rand = unifrnd(5) * unifrnd(3) + unifrnd(4)
233 : if (rand <= 0.0) cycle
234 : unifrnd(6) = (1._RKC + unifrnd(3)**2) * exp(unifrnd(4) * log(rand / unifrnd(4)) - unifrnd(5) * unifrnd(3))
235 : call setUnifRand(RNG unifrnd(7)) !call random number(unifrnd(7))
236 : if (unifrnd(7) <= unifrnd(6)) exit
237 : end do
238 : end if
239 : #if 0
240 : ! Alternative old implementation.
241 : function getRandGamma(alpha) result(randGamma)
242 : use pm_distNorm, only: setNormRand
243 : implicit none
244 : real(RK), intent(in) :: alpha
245 : real(RK) :: randGamma
246 : real(RK) :: c,u,v,z
247 : if (alpha<=0._RK) then ! illegal value of alpha
248 : randGamma = -1._RK
249 : return
250 : else
251 : randGamma = alpha
252 : if (randGamma < 1._RK) randGamma = randGamma + 1._RK
253 : randGamma = randGamma - 0.3333333333333333_RK
254 : c = 3._RK*sqrt(randGamma)
255 : c = 1._RK / c
256 : do
257 : do
258 : call setNormRand(z)
259 : v = 1._RK + c*z
260 : if (v<=0._RK) cycle
261 : exit
262 : end do
263 : v = v**3
264 : call setUnifRand(RNG u)
265 : if (log(u) >= 0.5_RK * z**2 + randGamma * (1._RK - v + log(v))) cycle
266 : randGamma = randGamma * v
267 : exit
268 : end do
269 : if (alpha < 1._RK) then
270 : call setUnifRand(RNG u)
271 : randGamma = randGamma * u**(1._RK/alpha)
272 : end if
273 : end if
274 : end function getRandGamma
275 : #endif
276 :
277 : #else
278 : !%%%%%%%%%%%%%%%%%%%%%%%%
279 : #error "Unrecognized interface."
280 : !%%%%%%%%%%%%%%%%%%%%%%%%
281 : #endif
282 : #undef GET_PARAM
283 : #undef GET_RAND
284 : #undef RNG
|
