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_distKolm](@ref pm_distKolm).
19 : !>
20 : !> \author
21 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
22 :
23 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 : #if getKolmPDF_ENABLED || setKolmPDF_ENABLED
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 :
29 : integer(IK) :: niter, niterSq
30 : real(TKC) :: expsq, temp, tempold, delta, xinv, exponent
31 : real(TKC), parameter :: xbreak = 1._TKC ! this seems to be the performance tipping point between the two series convergence in real64 and real32 precision.
32 : real(TKC), parameter :: SQRT2 = sqrt(2._TKC)
33 : real(TKC), parameter :: TWOSQRT2PI = 2 * SQRT2 * sqrt(acos(-1._TKC))
34 : real(TKC), parameter :: PI2SQRT8 = acos(-1._TKC) / sqrt(8._TKC)
35 : real(TKC), parameter :: TOL = 10 * epsilon(0._TKC)
36 : niter = 1_IK
37 : ! \devnote
38 : ! Based on experimentations, it takes roughly 2-4 term of the series in each branch to compute the PDF in [0, 3].
39 2020 : if (xbreak < x) then
40 1340 : exponent = SQRT2 * x
41 1340 : expsq = -exponent**2
42 1340 : pdf = exp(expsq) ! first term in the series.
43 1340 : if (pdf + TOL < 1._TKC) then
44 1336 : tempold = pdf
45 1336 : delta = 1._TKC
46 3960 : do
47 5302 : niter = niter + 1_IK
48 5302 : niterSq = niter * niter
49 5302 : temp = niterSq * exp(expsq * niterSq)
50 5302 : delta = sign(temp, -delta)
51 5302 : pdf = pdf + delta
52 5302 : if (tempold - temp < TOL) exit
53 3964 : tempold = temp
54 : end do
55 1340 : pdf = 8 * x * pdf
56 : else
57 0 : pdf = 0._TKC
58 : end if
59 : !print *, (niter + 1_IK) / 2_IK
60 680 : elseif (0._TKC < x) then ! warning: `x < 1` must hold.
61 674 : xinv = 1._TKC / x
62 674 : exponent = PI2SQRT8 * xinv
63 674 : expsq = exponent**2
64 674 : pdf = (expsq - .5_TKC) * exp(-expsq)
65 : do
66 1438 : niter = niter + 2_IK
67 1438 : expsq = (exponent * niter)**2
68 1438 : delta = (expsq - .5_TKC) * exp(-expsq)
69 1438 : pdf = pdf + delta
70 1438 : if (delta < TOL) exit
71 : end do
72 674 : pdf = pdf * TWOSQRT2PI * xinv**2
73 : !print *, (niter + 1_IK) / 2_IK
74 : else
75 6 : CHECK_ASSERTION(__LINE__, 0._TKC <= x, SK_"@setKolmPDF(): The condition `0 < x` must hold. x = "//getStr(x))
76 5 : pdf = 0._TKC
77 : end if
78 :
79 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 : #elif getKolmCDF_ENABLED || setKolmCDF_ENABLED
81 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 :
83 : integer(IK) :: niter
84 : real(TKC) :: temp, tempold, delta, xinv, exponent
85 : real(TKC), parameter :: xbreak = 1._TKC ! this seems to be the performance tipping point between the two series convergence in real64 and real32 precision.
86 : real(TKC), parameter :: SQRT2 = sqrt(2._TKC)
87 : real(TKC), parameter :: SQRT2PI = SQRT2 * sqrt(acos(-1._TKC))
88 : real(TKC), parameter :: PI2SQRT8 = acos(-1._TKC) / sqrt(8._TKC)
89 : real(TKC), parameter :: TOL = 10 * epsilon(0._TKC)
90 : niter = 1_IK
91 : ! \devnote
92 : ! Based on experimentations, it takes roughly two term of the series in each branch to compute the CDF.
93 2276 : if (xbreak < x) then
94 1447 : exponent = SQRT2 * x
95 1447 : cdf = exp(-exponent**2) ! first term in the series.
96 1447 : if (cdf + TOL < 1._TKC) then
97 1338 : tempold = cdf
98 1338 : delta = 1._TKC
99 3892 : do
100 5489 : niter = niter + 1_IK
101 5489 : temp = exp(-(exponent * niter)**2)
102 5489 : delta = sign(temp, -delta)
103 5489 : cdf = cdf + delta
104 5489 : if (tempold - temp < TOL) exit
105 4001 : tempold = temp
106 : end do
107 1447 : cdf = 1._TKC - 2 * cdf
108 : else
109 0 : cdf = 0._TKC
110 : end if
111 : !print *, (niter + 1_IK) / 2_IK
112 829 : elseif (0._TKC < x) then
113 819 : xinv = 1._TKC / x
114 819 : exponent = PI2SQRT8 * xinv
115 819 : cdf = exp(-exponent**2)
116 : do
117 1579 : niter = niter + 2_IK
118 1579 : delta = exp(-(exponent * niter)**2)
119 1579 : cdf = cdf + delta
120 1579 : if (delta < TOL) exit
121 : end do
122 819 : cdf = cdf * SQRT2PI * xinv
123 : !print *, (niter + 1_IK) / 2_IK
124 : else
125 10 : CHECK_ASSERTION(__LINE__, 0._TKC <= x, SK_"@setKolmPDF(): The condition `0 < x` must hold. x = "//getStr(x))
126 7 : cdf = 0._TKC
127 : end if
128 :
129 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 : #elif getKolmQuan_ENABLED || setKolmQuan_ENABLED
131 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 :
133 : real(TKC) :: cdfc
134 : real(TKC) :: logQuan, quanOld, func, ff, u, t
135 : real(TKC), parameter :: TOL = 100 * epsilon(0._TKC)
136 : real(TKC), parameter :: HALFPI = acos(-1._TKC) / 2._TKC
137 : real(TKC), parameter :: PIOVER8 = HALFPI / 4._TKC
138 6024 : CHECK_ASSERTION(__LINE__, 0._TKC <= cdf .and. cdf < 1._TKC, SK_"@setKolmQuan(): The condition `0 <= cdf .and. cdf < 1` must hold. cdf = "//getStr(cdf))
139 6024 : cdfc = 1._TKC - cdf
140 6024 : if (cdfc < 0.3_TKC) then
141 1843 : quan = 0.03_TKC
142 : do
143 23275 : quanOld = quan
144 23282 : quan = 0.5_TKC * cdfc + quan**4 - quan**9
145 23282 : if (quan > 0.06_TKC) quan = quan + quan**16 - quan**25
146 23282 : if (abs((quanOld - quan) / quan) < TOL) exit
147 : end do
148 1845 : quan = sqrt(-0.5_TKC * log(quan))
149 4179 : elseif (cdfc < 1._TKC) then
150 4169 : func = -PIOVER8 * (1._TKC - cdfc)**2
151 4169 : quan = getInvXlogX(func)
152 : do
153 21879 : logQuan = log(quan)
154 21879 : ff = func / (1._TKC + quan**4 + quan**12)**2
155 21879 : u = (quan * logQuan - ff) / (1._TKC + logQuan)
156 21879 : t = u / max(0.5_TKC, 1._TKC - 0.5_TKC * u / (quan * (1._TKC + logQuan)))
157 21879 : quan = quan - t
158 21879 : if (abs(t / quan) < TOL) exit
159 : end do
160 4169 : quan = HALFPI / sqrt(-log(quan))
161 : else
162 8 : quan = 0._TKC
163 : end if
164 :
165 : contains
166 :
167 4169 : PURE function getInvXlogX(y) result(invXlogX)
168 : real(TKC), parameter :: TOL = 10 * epsilon(0._TKC)
169 : real(TKC), parameter :: SQRTOL = sqrt(TOL)
170 : real(TKC), parameter :: STSTOL = sqrt(SQRTOL)
171 : real(TKC), parameter :: invNeper = exp(-1._TKC) ! 0.367879441171442322
172 : real(TKC), intent(in) :: y
173 : real(TKC) :: invXlogX
174 : real(TKC) :: t, to
175 4169 : CHECK_ASSERTION(__LINE__, -invNeper < y .and. y < 0._TKC, SK_"@setKolmQuan(): The condition `-exp(-1) < y .and. y < 0` must hold. y = "//getStr(y))
176 4163 : to = 0._TKC
177 4169 : if (y < -0.2) then
178 0 : invXlogX = log(invNeper - sqrt(2 * invNeper * (y + invNeper)))
179 : else
180 4163 : invXlogX = -10._TKC
181 : end if
182 22629 : do
183 26815 : t = (log(y / invXlogX) - invXlogX) * (invXlogX / (1._TKC + invXlogX))
184 26815 : invXlogX = invXlogX + t
185 26815 : if (t < SQRTOL .and. abs(t + to) < STSTOL * abs(t)) exit
186 26815 : if (abs(t / invXlogX) <= TOL) exit
187 26798 : to = t
188 : end do
189 4169 : invXlogX = exp(invXlogX)
190 4169 : end function
191 :
192 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
193 : #elif getKolmRand_ENABLED || setKolmRand_ENABLED
194 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195 :
196 3008 : CHECK_ASSERTION(__LINE__, 0._TKC <= unif .and. unif < 1._TKC, SK_"@setKolmQuan(): The condition `0 <= unif .and. unif < 1` must hold. unif = "//getStr(unif))
197 3008 : call setKolmQuan(rand, unif)
198 : #else
199 : !%%%%%%%%%%%%%%%%%%%%%%%%
200 : #error "Unrecognized interface."
201 : !%%%%%%%%%%%%%%%%%%%%%%%%
202 : #endif
|
