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 procedure implementations of [pm_distNorm](@ref pm_distNorm).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the RNG.
28 : #if setNormRand_ENABLED
29 : #if RNGD_ENABLED
30 : #define RNG
31 : #elif RNGF_ENABLED || RNGX_ENABLED
32 : #define RNG rng,
33 : #else
34 : #error "Unrecognized interface."
35 : #endif
36 : #endif
37 : !%%%%%%%%%%%%%%%%%%%%
38 : #if getNormLogPDF_ENABLED
39 : !%%%%%%%%%%%%%%%%%%%%
40 :
41 : real(RKC) :: mu_def, invSigma
42 14968 : if (present(mu)) then
43 4018 : mu_def = mu
44 : else
45 10950 : mu_def = 0._RKC
46 : end if
47 14968 : if (present(sigma)) then
48 3016 : invSigma = 1._RKC / sigma
49 : else
50 11952 : invSigma = 1._RKC
51 : end if
52 14968 : call setNormLogPDF(logPDF, x, mu_def, invSigma, log(invSigma))
53 :
54 : !%%%%%%%%%%%%%%%%%%%%
55 : #elif setNormLogPDF_ENABLED
56 : !%%%%%%%%%%%%%%%%%%%%
57 :
58 : real(RKC) , parameter :: PI = acos(-1._RKC)
59 : real(RKC) , parameter :: INVERSE_SQRT_TWO_PI = 1._RKC / sqrt(2._RKC * PI)
60 : real(RKC) , parameter :: LOG_INVERSE_SQRT_TWO_PI = log(INVERSE_SQRT_TWO_PI)
61 : #if DS_ENABLED || MS_ENABLED
62 34889 : CHECK_ASSERTION(__LINE__, 0._RKC < invSigma, SK_"@setNormLogPDF(): The condition `0. < invSigma` must hold. invSigma = "//getStr(invSigma))
63 104667 : CHECK_ASSERTION(__LINE__, abs(logInvSigma - log(invSigma)) <= epsilon(0._RKC)*100, SK_"@setNormLogPDF(): The condition `abs(logInvSigma - log(invSigma)) <= epsilon(0._RKC)*100` must hold. logInvSigma, log(invSigma) = "//getStr([logInvSigma, log(invSigma)]))
64 : #endif
65 : #if DD_ENABLED
66 78990 : logPDF = LOG_INVERSE_SQRT_TWO_PI - 0.5_RKC * x**2
67 : #elif DS_ENABLED
68 1 : logPDF = logInvSigma + LOG_INVERSE_SQRT_TWO_PI - 0.5_RKC * (x * invSigma)**2
69 : #elif MD_ENABLED
70 1 : logPDF = LOG_INVERSE_SQRT_TWO_PI - 0.5_RKC * (x - mu)**2
71 : #elif MS_ENABLED
72 34888 : logPDF = logInvSigma + LOG_INVERSE_SQRT_TWO_PI - 0.5_RKC * ((x - mu) * invSigma)**2
73 : #else
74 : #error "Unrecognized interface"
75 : #endif
76 :
77 : !%%%%%%%%%%%%%%%%%
78 : #elif getNormCDF_ENABLED
79 : !%%%%%%%%%%%%%%%%%
80 :
81 4522 : if (present(mu) .and. present(sigma)) then
82 4024 : call setNormCDF(cdf, x, mu, 1._RKC / sigma)
83 498 : elseif (present(sigma)) then
84 0 : call setNormCDF(cdf, x, 0._RKC, 1._RKC / sigma)
85 498 : elseif (present(mu)) then
86 2 : call setNormCDF(cdf, x, mu)
87 : else
88 496 : call setNormCDF(cdf, x)
89 : end if
90 :
91 : !%%%%%%%%%%%%%%%%%
92 : #elif setNormCDF_ENABLED
93 : !%%%%%%%%%%%%%%%%%
94 :
95 : real(RKC) , parameter :: INVERSE_SQRT_TWO = 1._RKC / sqrt(2._RKC)
96 : #if DD_ENABLED
97 840 : cdf = 0.5_RKC * (1._RKC + erf(x * INVERSE_SQRT_TWO))
98 : #elif MD_ENABLED
99 4 : cdf = 0.5_RKC * (1._RKC + erf((x - mu) * INVERSE_SQRT_TWO))
100 : #elif MS_ENABLED
101 8040 : CHECK_ASSERTION(__LINE__, 0._RKC < invSigma, SK_"@setNormCDF(): The condition `0. < invSigma` must hold. invSigma = "//getStr(invSigma))
102 8040 : cdf = 0.5_RKC * (1._RKC + erf((x - mu) * invSigma * INVERSE_SQRT_TWO))
103 : #else
104 : #error "Unrecognized interface"
105 : #endif
106 :
107 : !%%%%%%%%%%%%%%%%%%
108 : #elif getNormQuan_ENABLED
109 : !%%%%%%%%%%%%%%%%%%
110 :
111 4021 : if (present(mu) .and. present(sigma)) then
112 4007 : call setNormQuan(quantile, cdf, mu, sigma)
113 14 : elseif (present(sigma)) then
114 0 : call setNormQuan(quantile, cdf, 0._RKC, sigma)
115 14 : elseif (present(mu)) then
116 7 : call setNormQuan(quantile, cdf, mu)
117 : else
118 7 : call setNormQuan(quantile, cdf)
119 : end if
120 :
121 : !%%%%%%%%%%%%%%%%%%
122 : #elif setNormQuan_ENABLED
123 : !%%%%%%%%%%%%%%%%%%
124 :
125 : real(RKC), parameter :: SQRT_TWO = sqrt(2._RKC)
126 8042 : if (0._RKC < cdf .and. cdf < 1._RKC) then
127 8030 : call setErfInv(quantile, 2 * cdf - 1._RKC, abserr = 100 * epsilon(cdf))
128 : #if DD_ENABLED
129 10 : quantile = quantile * SQRT_TWO
130 : #elif MD_ENABLED
131 10 : quantile = quantile * SQRT_TWO + mu
132 : #elif MS_ENABLED
133 8010 : CHECK_ASSERTION(__LINE__, 0._RKC < sigma, SK_"@setNormQuan(): The condition `0. < sigma` must hold. sigma = "//getStr(sigma))
134 8010 : quantile = quantile * SQRT_TWO * sigma + mu
135 : #else
136 : #error "Unrecognized interface"
137 : #endif
138 : else
139 12 : quantile = huge(quantile)
140 12 : if (0._RKC == cdf) quantile = -quantile
141 12 : CHECK_ASSERTION(__LINE__, 0._RKC <= cdf .and. cdf <= 1._RKC, SK_"@setNormQuan(): The condition `0. <= cdf .and. cdf <= 1` must hold. cdf = "//getStr(cdf))
142 : end if
143 :
144 : !%%%%%%%%%%%%%%%%%
145 : #elif getZigNorm_ENABLED
146 : !%%%%%%%%%%%%%%%%%
147 :
148 : zig = getZig( nlay = nlay & ! LCOV_EXCL_LINE
149 : , getFunc = getFuncNorm_RKC & ! LCOV_EXCL_LINE
150 : , getFuncInv = getFuncInvNorm_RKC & ! LCOV_EXCL_LINE
151 : , getZigArea = getZigAreaNorm_RKC & ! LCOV_EXCL_LINE
152 : , abserr = abserr & ! LCOV_EXCL_LINE
153 : , abstol = abstol & ! LCOV_EXCL_LINE
154 3 : )
155 :
156 : contains
157 :
158 7796 : PURE function getFuncNorm_RKC(x) result(func)
159 : real(RKC) , intent(in) :: x
160 : real(RKC) :: func
161 : !func = getFuncNorm(x)
162 7862 : call setNormLogPDF(func, x)
163 7862 : func = exp(func)
164 7796 : end function
165 :
166 7173 : pure function getFuncInvNorm_RKC(func) result(x)
167 : real(RKC) , parameter :: LOGSQRT2PI = log(sqrt(2 * acos(-1._RKC)))
168 : real(RKC) , intent(in) :: func
169 : real(RKC) :: x
170 7173 : x = sqrt(-2._RKC * (LOGSQRT2PI + log(func)))
171 : !x = getFuncInvNorm(func)
172 7173 : end function
173 :
174 343 : PURE function getZigAreaNorm_RKC(r) result(area)
175 : real(RKC) , intent(in) :: r
176 : real(RKC) :: area
177 343 : CHECK_ASSERTION(__LINE__, 0._RKC <= r, SK_"@getZigAreaNorm(): The condition `0. <= r` must hold. r = "//getStr(r))
178 : !area = getZigAreaNorm(r)
179 343 : call setNormCDF(area, r)
180 343 : area = 1._RKC - area + r * getFuncNorm_RKC(r)
181 343 : end function
182 :
183 : !pure function getGradNorm_RKC(x) result(grad)
184 : ! real(RKC) , intent(in) :: x
185 : ! real(RKC) :: grad
186 : ! grad = getGradNorm(x)
187 : !end function
188 :
189 : ! !%%%%%%%%%%%%%%%%%%
190 : !#elif getFuncNorm_ENABLED
191 : ! !%%%%%%%%%%%%%%%%%%
192 : !
193 : ! real(RKC), parameter :: INVSQRT2PI = 1._RKC / sqrt(2 * acos(-1._RKC))
194 : ! func = INVSQRT2PI * exp(-0.5_RKC * x**2)
195 : !
196 : ! !%%%%%%%%%%%%%%%%%%
197 : !#elif getGradNorm_ENABLED
198 : ! !%%%%%%%%%%%%%%%%%%
199 : !
200 : ! grad = -x * getFuncNorm(x)
201 : !
202 : ! !%%%%%%%%%%%%%%%%%%%%%
203 : !#elif getFuncInvNorm_ENABLED
204 : ! !%%%%%%%%%%%%%%%%%%%%%
205 : !
206 : ! x = sqrt(-2._RKC * log(func))
207 : !
208 : ! !%%%%%%%%%%%%%%%%%%%%%
209 : !#elif getZigAreaNorm_ENABLED
210 : ! !%%%%%%%%%%%%%%%%%%%%%
211 : !
212 : ! real(RKC), parameter :: SQRT2PI = sqrt(2 * acos(-1._RKC))
213 : ! CHECK_ASSERTION(__LINE__, 0._RKC <= r, SK_"@getZigAreaNorm(): The condition `0. <= r` must hold. r = "//getStr(r))
214 : ! call setNormCDF(area, r)
215 : ! area = r * getFuncNorm(r) + SQRT2PI * (1._RKC - area)
216 :
217 : !%%%%%%%%%%%%%%%%%%
218 : #elif getNormRand_ENABLED
219 : !%%%%%%%%%%%%%%%%%%
220 :
221 15168 : call setNormRand(rand)
222 15168 : if (present(std)) rand = rand * std
223 15168 : rand = rand + mean
224 :
225 : !%%%%%%%%%%%%%%%%%%
226 : #elif setNormRand_ENABLED
227 : !%%%%%%%%%%%%%%%%%%
228 :
229 : integer(IK) :: index
230 : real(RKC) :: dumm, func
231 : #if ZD_ENABLED
232 : real(RKC) , parameter :: zig(1 : 2, 0 : 256) = real(ZIG_RKB, RKC)
233 : integer(IK) , parameter :: nlay = ubound(zig, 2, IK)
234 3916538 : CHECK_ASSERTION(__LINE__, precision(rand) <= ZIG_PRECISION, SK_"@setNormRand(): The condition `precision(rand) <= ZIG_PRECISION` must hold. precision(rand) <= ZIG_PRECISION = "//getStr([int(precision(rand), IK), ZIG_PRECISION]))
235 : #elif ZA_ENABLED
236 : integer(IK) :: nlay
237 2 : nlay = ubound(zig, 2, IK)
238 20 : CHECK_ASSERTION(__LINE__, all(0._RKC <= zig), SK_"@setNormRand(): The condition `all(0. <= zig)` must hold. zig = "//getStr(zig))
239 2 : CHECK_ASSERTION(__LINE__, size(zig, 1, IK) == 2_IK, SK_"@setNormRand(): The condition `size(zig, 1) == 2` must hold. shape(zig) = "//getStr(shape(zig, IK)))
240 : #else
241 : #error "Unrecognized interface."
242 : #endif
243 : #if D1_ENABLED
244 : #define GET_RAND(i)rand(i)
245 : block
246 : integer(IK) :: irand
247 86283 : loopOverRand: do irand = 1_IK, size(rand, 1, IK)
248 : #elif D0_ENABLED
249 : #define GET_RAND(i)rand
250 : #else
251 : #error "Unrecognized interface."
252 : #endif
253 15052 : loopTry: do
254 4004124 : call setUnifRand(RNG index, 0_IK, nlay - 1)
255 : #if RNGD_ENABLED
256 : #define SET_RAND(X)call random_number(X)
257 254371 : call random_number(GET_RAND(irand))
258 : !if (GET_RAND(irand) == 0 .or. GET_RAND(irand) == 1) print *, "GET_RAND(irand)", GET_RAND(irand)
259 254371 : GET_RAND(irand) = GET_RAND(irand) * 2 - 1._RKC
260 : #elif RNGF_ENABLED || RNGX_ENABLED
261 : #define SET_RAND(X)call setUnifRand(rng, X)
262 3749753 : call setUnifRand(RNG GET_RAND(irand), -1._RKC, +1._RKC)
263 : #else
264 : #error "Unrecognized interface."
265 : #endif
266 4004124 : GET_RAND(irand) = GET_RAND(irand) * zig(1, index)
267 4004124 : if (zig(1, index + 1) < abs(GET_RAND(irand))) then ! rejection involved.
268 75685 : SET_RAND(dumm)
269 75685 : if (0_IK < index) then
270 71127 : dumm = zig(2, index) + dumm * (zig(2, index + 1) - zig(2, index)) ! height
271 71127 : call setNormLogPDF(func, GET_RAND(irand)); func = exp(func)
272 71127 : if (func < dumm) cycle loopTry
273 : else
274 : ! For the normal tail, we follow the method of Marsaglia, 1964:
275 : ! generate x = -ln(U1) = r; y = -ln(U2) until y + y > x * x, then return r + x.
276 4558 : index = int(sign(1._RKC, GET_RAND(irand)), IK)
277 1331 : loopTail: do
278 : !print *, "loopTail", dumm, zig(1)
279 5889 : GET_RAND(irand) = -log(1._RKC - dumm) / zig(1, 1)
280 5889 : SET_RAND(dumm)
281 : !if (dumm == 0 .or. dumm == 1) print *, "dumm2", dumm
282 5889 : dumm = -log(1._RKC - dumm)
283 5889 : if (GET_RAND(irand)**2 < 2 * dumm) exit loopTail
284 1331 : SET_RAND(dumm)
285 : !if (dumm == 0 .or. dumm == 1) print *, "dumm3", dumm
286 : end do loopTail
287 4558 : GET_RAND(irand) = GET_RAND(irand) + zig(1, 1)
288 4558 : if (index < 0_IK) GET_RAND(irand) = -GET_RAND(irand)
289 : end if
290 : end if
291 31405 : exit loopTry
292 : end do loopTry
293 : #if D1_ENABLED
294 : end do loopOverRand
295 : end block
296 : #endif
297 : #undef SET_RAND
298 : #undef GET_RAND
299 :
300 : !%%%%%%%%%%%%%%%%%%%%%
301 : #elif getNormRandBox_ENABLED
302 : !%%%%%%%%%%%%%%%%%%%%%
303 :
304 : ! Box-Muller rejection approach.
305 : real(RKC) :: rand1
306 : real(RKC) , save :: rand2 = -1._RKC
307 : logical(LK) , save :: failed = .true._LK
308 : if (failed) then
309 : do
310 : call random_number(rand1)
311 : call random_number(rand2)
312 : call setNormRandBox(rand1, rand2, failed)
313 : if (.not. failed) exit
314 : end do
315 : rand = rand1
316 : else
317 : rand = rand2
318 : failed = .true._LK
319 : end if
320 : #if MD_ENABLED
321 : rand = rand + mean
322 : #elif MS_ENABLED
323 : rand = rand * std + mean
324 : #else
325 : #error "Unrecognized interface."
326 : #endif
327 :
328 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
329 : #elif setNormRandBox_ENABLED && D0_ENABLED
330 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
331 :
332 : ! Perform checks.
333 : #define CHECK_RAND1 \
334 : CHECK_ASSERTION(__LINE__, 0._RKC <= rand1 .and. rand1 < 1._RKC, SK_"@setNormRand(): The condition `0. <= rand1 .and. rand1 < 1.` must hold. rand1 = "//getStr(rand1));
335 : #define CHECK_RAND2 \
336 : CHECK_ASSERTION(__LINE__, 0._RKC <= rand2 .and. rand2 < 1._RKC, SK_"@setNormRand(): The condition `0. <= rand2 .and. rand2 < 1.` must hold. rand2 = "//getStr(rand2));
337 : ! Define the factor.
338 : #if Basic_ENABLED
339 : #define GET_RE(x) x%re
340 : #define GET_IM(x) x%im
341 : complex(RKC) :: factor
342 : real(RKC), parameter :: TWP_PI = 2 * acos(-1._RKC), EPS = epsilon(0._RKC)
343 7501 : CHECK_RAND1
344 7501 : CHECK_RAND2
345 7501 : factor = sqrt(-2._RKC * log(rand1 + EPS)) * exp(cmplx(0._RKC, -TWP_PI * rand2, RKC)) ! compute sin and cos, hopefully in one instruction.
346 : #elif Polar_ENABLED
347 : #define GET_RE(x) rand1 * x
348 : #define GET_IM(x) rand2 * x
349 : real(RKC) :: factor, rsq
350 19097 : CHECK_RAND1
351 19097 : CHECK_RAND2
352 19097 : rand1 = 2._RKC * rand1 - 1._RKC
353 19097 : rand2 = 2._RKC * rand2 - 1._RKC
354 19097 : rsq = rand1**2 + rand2**2
355 19097 : if (0._RKC < rsq .and. rsq < 1._RKC) then
356 15000 : failed = .false._LK
357 15000 : factor = sqrt(-2._RKC * log(rsq) / rsq)
358 : else
359 4097 : failed = .true._LK
360 4097 : return
361 : end if
362 : #else
363 : #error "Unrecognized interface."
364 : #endif
365 : ! Scale the numbers.
366 : #if MS_ENABLED
367 : CHECK_ASSERTION(__LINE__, 0._RKC < std, SK_"@setNormRand(): The condition `0. < std` must hold. std = "//getStr(std))
368 : rand1 = GET_RE(factor) * std + mean
369 : rand2 = GET_IM(factor) * std + mean
370 : #elif MD_ENABLED
371 : rand1 = GET_RE(factor) + mean
372 : rand2 = GET_IM(factor) + mean
373 : #elif DD_ENABLED
374 22501 : rand1 = GET_RE(factor)
375 15000 : rand2 = GET_IM(factor)
376 : #else
377 : #error "Unrecognized interface."
378 : #endif
379 :
380 : !%%%%%%%%%%%%%%%%%%%%%
381 : #elif getNormEntropy_ENABLED
382 : !%%%%%%%%%%%%%%%%%%%%%
383 :
384 : real(RKC), parameter :: LOG_TWO_PI = log(2 * acos(-1._RKC))
385 7 : entropy = 0.5_RKC + 0.5_RKC * (LOG_TWO_PI + logVar)
386 :
387 : !%%%%%%%%%%%%%%%%%%%%
388 : #elif getNormFisher_ENABLED
389 : !%%%%%%%%%%%%%%%%%%%%
390 :
391 4 : CHECK_ASSERTION(__LINE__, 0._RKC <= varinv, SK_"@getNormFisher(): The input `varInv` must be positive. varInv = "//getStr(varInv))
392 4 : Fisher(1,1) = varInv
393 4 : Fisher(2,1) = 0._RKC
394 4 : Fisher(1,2) = 0._RKC
395 4 : Fisher(2,2) = 2._RKC * varInv
396 :
397 : !%%%%%%%%%%%%%%%%%
398 : #elif getNormKLD_ENABLED
399 : !%%%%%%%%%%%%%%%%%
400 :
401 : #if DV_ENABLED || MV_ENABLED
402 : real(RKC) :: varP2varQ
403 9 : CHECK_ASSERTION(__LINE__, 0._RKC < varP, SK_"@getNormKLD(): The condition `0. < varP` must hold. varP = "//getStr(varP))
404 9 : CHECK_ASSERTION(__LINE__, 0._RKC < varQ, SK_"@getNormKLD(): The condition `0. < varQ` must hold. varP = "//getStr(varQ))
405 : #endif
406 : #if MD_ENABLED || MV_ENABLED
407 6 : CHECK_ASSERTION(__LINE__, 0._RKC <= meanDiffSq, SK_"@getNormKLD(): The condition `0. <= meanDiffSq` must hold. meanDiffSq = "//getStr(meanDiffSq))
408 : #endif
409 : #if MV_ENABLED
410 3 : varP2varQ = varP / varQ
411 3 : kld = 0.5_RKC * (varP2varQ - log(varP2varQ) - 1._RKC + meanDiffSq / varQ)
412 : #elif MD_ENABLED
413 3 : kld = 0.5_RKC * meanDiffSq
414 : #elif DV_ENABLED
415 6 : varP2varQ = varP / varQ
416 6 : kld = 0.5_RKC * (varP2varQ - log(varP2varQ) - 1._RKC)
417 : #else
418 : #error "Unrecognized interface."
419 : #endif
420 :
421 : #else
422 : !%%%%%%%%%%%%%%%%%%%%%%%%
423 : #error "Unrecognized interface."
424 : !%%%%%%%%%%%%%%%%%%%%%%%%
425 : #endif
426 :
427 : #undef CHECK_RAND1
428 : #undef CHECK_RAND2
429 : #undef GET_RE
430 : #undef GET_IM
431 : #undef RNG
|
