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_distPiwiPoweto](@ref pm_distPiwiPoweto).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getPiwiPowetoLogPDFNF_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : use pm_mathCumSum, only: setCumSum
32 : use pm_mathLogSubExp, only: getLogSubExp
33 : use pm_mathLogSumExp, only: getLogSumExp
34 : real(RKC), parameter :: LOG_HUGE = log(huge(logLimX))
35 : integer(IK) :: i, lenAlpha, lenLogLimX
36 : real(RKC) :: maxArea, logIntegral
37 : #if ALD_ENABLED
38 7466 : real(RKC) :: cumSumArea(size(logLimX, kind = IK)) ! initially, LogArea.
39 : #elif ALC_ENABLED
40 53871 : CHECK_ASSERTION(__LINE__, size(logLimX,1,IK) == size(cumSumArea,1,IK), SK_"@getPiwiPowetoLogPDFNF(): The condition `size(logLimX,1,IK) == size(cumSumArea,1,IK)` must hold. size(logLimX), size(cumSumArea) = "//getStr([size(logLimX,1,IK), size(cumSumArea,1,IK)]))
41 : #else
42 : #error "Unrecognized interface."
43 : #endif
44 25423 : lenAlpha = size(alpha, kind = IK)
45 : lenLogLimX = size(logLimX, kind = IK)
46 25423 : CHECK_ASSERTION(__LINE__, lenAlpha > 0_IK, SK_"@getPiwiPowetoLogPDFNF(): The condition `size(alpha) > 0` must hold. size(alpha) = "//getStr(lenAlpha))
47 76269 : CHECK_ASSERTION(__LINE__, lenLogLimX == lenAlpha + 1_IK, SK_"@getPiwiPowetoLogPDFNF(): The condition `size(logLimX) == size(alpha) + 1` must hold. size(logLimX), size(alpha) = "//getStr([lenLogLimX, lenAlpha]))
48 25423 : CHECK_ASSERTION(__LINE__, isAscending(logLimX), SK_"@getPiwiPowetoLogPDFNF(): The condition `isAscending(logLimX)` must hold. logLimX = "//getStr(logLimX))
49 101692 : CHECK_ASSERTION(__LINE__, alpha(1) > 0._RKC .or. logLimX(1) > -Log_HUGE, SK_"@getPiwiPowetoLogPDFNF(): The conditions `alpha(1) > 0._RKC .or. logLimX(1) > -Log_HUGE` must hold. alpha(1), logLimX(1), -Log_HUGE = "//getStr([alpha(1), logLimX(1), -Log_HUGE]))
50 101692 : CHECK_ASSERTION(__LINE__, alpha(lenAlpha) < 0._RKC .or. logLimX(lenLogLimX) < Log_HUGE, SK_"@getPiwiPowetoLogPDFNF(): The conditions `alpha(lenAlpha) < 0._RKC .or. logLimX(lenLogLimX) < Log_HUGE` must hold. alpha(lenAlpha), logLimX(lenLogLimX), Log_HUGE = "//getStr([alpha(lenAlpha), logLimX(lenLogLimX), Log_HUGE]))
51 :
52 : ! Initially, cumSumArea contains areas of the segments.
53 25423 : if (alpha(lenAlpha) > 0._RKC) then ! Power tail.
54 146 : cumSumArea(lenLogLimX) = getLogSubExp(smaller = alpha(lenAlpha) * logLimX(lenAlpha), larger = alpha(lenAlpha) * logLimX(lenLogLimX)) - log(alpha(lenAlpha))
55 25277 : elseif (alpha(lenAlpha) < 0._RKC) then
56 23178 : cumSumArea(lenLogLimX) = getLogSubExp(smaller = alpha(lenAlpha) * logLimX(lenLogLimX), larger = alpha(lenAlpha) * logLimX(lenAlpha)) - log(-alpha(lenAlpha))
57 : else
58 2099 : cumSumArea(lenLogLimX) = log(logLimX(lenLogLimX) - logLimX(lenAlpha))
59 : end if
60 25423 : logPDFNF(lenAlpha) = 0._RKC
61 25423 : maxArea = cumSumArea(lenLogLimX)
62 92753 : do i = lenAlpha - 1, 1_IK, -1_IK
63 67330 : logPDFNF(i) = logPDFNF(i + 1) + (alpha(i + 1) - alpha(i)) * logLimX(i + 1)
64 67330 : if (alpha(i) > 0._RKC) then
65 33393 : cumSumArea(i + 1) = logPDFNF(i) + getLogSubExp(smaller = alpha(i) * logLimX(i), larger = alpha(i) * logLimX(i + 1)) - log(alpha(i))
66 33937 : elseif (alpha(i) < 0._RKC) then
67 25251 : cumSumArea(i + 1) = logPDFNF(i) + getLogSubExp(smaller = alpha(i) * logLimX(i + 1), larger = alpha(i) * logLimX(i)) - log(-alpha(i))
68 : else
69 8686 : cumSumArea(i + 1) = logPDFNF(i) + log(logLimX(i + 1) - logLimX(i))
70 : end if
71 92753 : if (maxArea < cumSumArea(i + 1)) maxArea = cumSumArea(i + 1)
72 : end do
73 25423 : logIntegral = getLogSumExp(cumSumArea(2:lenLogLimX), maxArea)
74 118176 : logPDFNF = logPDFNF - logIntegral
75 : #if ALC_ENABLED
76 17957 : cumSumArea(1) = 0._RKC
77 82241 : cumSumArea(2:lenLogLimX) = exp(cumSumArea(2:lenLogLimX) - logIntegral)
78 17957 : call setCumSum(cumSumArea(2:lenLogLimX))
79 : #endif
80 :
81 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
82 : #elif getPiwiPowetoLogPDF_ENABLED
83 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
84 :
85 : if (present(logPDFNF)) then
86 81 : call setPiwiPowetoLogPDF(logPDF, logx, alpha, logLimX, logPDFNF)
87 : else
88 2487 : call setPiwiPowetoLogPDF(logPDF, logx, alpha, logLimX, getPiwiPowetoLogPDFNF(alpha, logLimX))
89 : end if
90 :
91 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92 : #elif setPiwiPowetoLogPDF_ENABLED && ALL_ENABLED
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : integer(IK) :: i, lenAlpha
96 92001 : lenAlpha = size(alpha, kind = IK)
97 92001 : CHECK_ASSERTION(__LINE__, lenAlpha > 0_IK, SK_"@setPiwiPowetoLogPDF(): The condition `size(alpha) > 0` must hold. size(logLimX) = "//getStr(lenAlpha))
98 276003 : CHECK_ASSERTION(__LINE__, size(logPDFNF,1,IK) == lenAlpha, SK_"@setPiwiPowetoLogPDF(): The condition `size(logPDFNF,1,IK) == lenAlpha` must hold. size(logPDFNF,1,IK), lenAlpha = "//getStr([size(logPDFNF,1,IK), lenAlpha]))
99 276003 : CHECK_ASSERTION(__LINE__, size(logLimX,1,IK) == lenAlpha + 1_IK, SK_"@setPiwiPowetoLogPDF(): The condition `size(logLimX) == size(alpha) + 1` must hold. size(logLimX,1,IK), size(alpha,1,IK) = "//getStr([size(logLimX,1,IK), size(alpha,1,IK)]))
100 92001 : CHECK_ASSERTION(__LINE__, isAscending(logLimX), SK_"@setPiwiPowetoLogPDF(): The conditions `isAscending(logLimX)` must hold. logLimX = "//getStr(logLimX))
101 276003 : CHECK_ASSERTION(__LINE__, logLimX(1) <= logx, SK_"@setPiwiPowetoLogPDFNF(): The condition `logLimX(1) <= logx` must hold. logLimX(1), logx = "//getStr([logLimX(1), logx]))
102 276003 : CHECK_ASSERTION(__LINE__, logx <= logLimX(size(logLimX,1,IK)), SK_"@getPiwiPowetoLogPDFNF(): The condition `logx < logLimX(size(logLimX,1,IK))` must hold. logx, logLimX(size(logLimX,1,IK)) = "//getStr([logx, logLimX(size(logLimX,1,IK))]))
103 : !CHECK_ASSERTION(__LINE__, logLimX(size(logLimX,1,IK)) < log(huge(logLimX)), SK_"@getPiwiPowetoLogPDFNF(): The condition `logLimX(size(logLimX)) < log(huge(logLimX))` must hold. logLimX = "//getStr(logLimX))
104 182600 : do i = lenAlpha, 1_IK, -1_IK
105 182600 : if (logx < logLimX(i)) cycle
106 92001 : logPDF = logPDFNF(i) + (alpha(i) - 1._RKC) * logx
107 182600 : return
108 : end do
109 :
110 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111 : #elif setPiwiPowetoLogPDF_ENABLED && BAN_ENABLED
112 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113 :
114 243 : CHECK_ASSERTION(__LINE__, 0_IK < bin .and. bin <= size(alpha, kind = IK), SK_"@setPiwiPowetoLogPDFNF(): The conditions `0_IK < bin .and. bin <= size(alpha, kind = IK)` must hold. bin, size(alpha, kind = IK) = "//getStr([bin, size(alpha, kind = IK)]))
115 243 : CHECK_ASSERTION(__LINE__, size(alpha, kind = IK) == size(logPDFNF, kind = IK), SK_"@setPiwiPowetoLogPDFNF(): The condition `size(alpha) == size(logPDFNF)` must hold. size(alpha), size(logPDFNF) = "//getStr([size(alpha, kind = IK), size(logPDFNF, kind = IK)]))
116 81 : logPDF = logPDFNF(bin) + (alpha(bin) - 1._RKC) * logx
117 :
118 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 : #elif getPiwiPowetoCDF_ENABLED && ALDD_ENABLED
120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 :
122 1478 : real(RKC) :: cumSumArea(size(logLimX, kind = IK))
123 739 : call setPiwiPowetoCDF(cdf, logx, alpha, logLimX, getPiwiPowetoLogPDFNF(alpha, logLimX, cumSumArea), cumSumArea)
124 :
125 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126 : #elif getPiwiPowetoCDF_ENABLED && ALLC_ENABLED
127 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 :
129 2805 : call setPiwiPowetoCDF(cdf, logx, alpha, logLimX, logPDFNF, cumSumArea)
130 :
131 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132 : #elif setPiwiPowetoCDF_ENABLED && BAN_ENABLED
133 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134 :
135 : use pm_mathLogSubExp, only: getLogSubExp
136 34812 : CHECK_ASSERTION(__LINE__, size(logLimX, kind = IK) == size(cumSumArea, kind = IK), SK_"@setPiwiPowetoCDF(): The condition `size(logLimX) == size(cumSumArea)` must hold. size(logLimX), size(cumSumArea) = "//getStr([size(logLimX, kind = IK), size(cumSumArea, kind = IK)]))
137 34812 : CHECK_ASSERTION(__LINE__, size(logLimX, kind = IK) == size(alpha, kind = IK) + 1_IK, SK_"@setPiwiPowetoCDF(): The condition `size(logLimX) == size(alpha) + 1` must hold. size(logLimX), size(alpha) = "//getStr([size(logLimX, kind = IK), size(alpha, kind = IK)]))
138 34812 : CHECK_ASSERTION(__LINE__, size(alpha, kind = IK) == size(logPDFNF, kind = IK), SK_"@setPiwiPowetoCDF(): The condition `size(alpha) == size(logPDFNF)` must hold. size(alpha), size(logPDFNF) = "//getStr([size(alpha, kind = IK), size(logPDFNF, kind = IK)]))
139 34812 : CHECK_ASSERTION(__LINE__, 0_IK < bin .and. bin < size(logLimX, kind = IK), SK_"@setPiwiPowetoCDF(): The conditions `0_IK < bin .and. bin < size(logLimX)` must hold. bin, size(logLimX) = "//getStr([bin, size(logLimX, kind = IK)]))
140 46416 : CHECK_ASSERTION(__LINE__, logLimX(bin) <= logx .and. logx <= logLimX(bin + 1), SK_"@setPiwiPowetoCDF(): The conditions `logLimX(bin) <= logx .and. logx < logLimX(bin + 1)` must hold. logLimX(bin), logx, logLimX(bin+1) = "//getStr([logLimX(bin), logx, logLimX(bin+1)]))
141 : #if CHECK_ENABLED
142 : block
143 23208 : real(RKC) :: logPDFNF_ref(size(logPDFNF, kind = IK)), cumSumArea_ref(size(cumSumArea, kind = IK))
144 11604 : logPDFNF_ref = getPiwiPowetoLogPDFNF(alpha, logLimX, cumSumArea_ref)
145 134277 : CHECK_ASSERTION(__LINE__, all(abs(logPDFNF - logPDFNF_ref) <= epsilon(0._RKC) * 1000), SK_"@setPiwiPowetoCDF(): The condition `all(abs(logPDFNF - logPDFNF_ref) <= epsilon(0._RKC) * 1000)` must hold. logPDFNF - logPDFNF_ref = "//getStr([logPDFNF - logPDFNF_ref]))
146 169089 : CHECK_ASSERTION(__LINE__, all(abs(cumSumArea - cumSumArea_ref) <= epsilon(0._RKC) * 1000), SK_"@setPiwiPowetoCDF(): The condition `all(abs(cumSumArea - cumSumArea_ref) <= epsilon(0._RKC) * 1000)` must hold. cumSumArea - cumSumArea_ref = "//getStr([cumSumArea - cumSumArea_ref]))
147 : end block
148 : #endif
149 11604 : cdf = cumSumArea(bin)
150 11604 : if (logx /= logLimX(bin)) then
151 8814 : if (alpha(bin) > 0._RKC) then
152 4593 : cdf = cumSumArea(bin) + exp(logPDFNF(bin) + getLogSubExp(smaller = alpha(bin) * logLimX(bin), larger = alpha(bin) * logx)) / alpha(bin)
153 4221 : elseif (alpha(bin) < 0._RKC) then
154 3325 : cdf = cumSumArea(bin) - exp(logPDFNF(bin) + getLogSubExp(smaller = alpha(bin) * logx, larger = alpha(bin) * logLimX(bin))) / alpha(bin)
155 : else
156 896 : cdf = cumSumArea(bin) + exp(logPDFNF(bin) + log(logx - logLimX(bin)))
157 : end if
158 : end if
159 :
160 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161 : #elif setPiwiPowetoCDF_ENABLED && MAN_ENABLED
162 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163 :
164 : integer(IK) :: i, lenAlpha
165 9176 : lenAlpha = size(alpha, kind = IK)
166 9176 : CHECK_ASSERTION(__LINE__, lenAlpha > 0_IK, SK_"@setPiwiPowetoCDF(): The condition `size(alpha) > 0` must hold. size(logLimX) = "//getStr(lenAlpha))
167 27528 : CHECK_ASSERTION(__LINE__, size(logPDFNF,1,IK) == lenAlpha, SK_"@setPiwiPowetoCDF(): The condition `size(logPDFNF,1,IK) == lenAlpha` must hold. size(logPDFNF,1,IK), lenAlpha = "//getStr([size(logPDFNF,1,IK), lenAlpha]))
168 27528 : CHECK_ASSERTION(__LINE__, size(logLimX,1,IK) == lenAlpha + 1_IK, SK_"@setPiwiPowetoCDF(): The condition `size(logLimX) == size(alpha) + 1` must hold. size(logLimX,1,IK), size(alpha,1,IK) = "//getStr([size(logLimX,1,IK), size(alpha,1,IK)]))
169 36704 : CHECK_ASSERTION(__LINE__, logLimX(1) <= logx .and. logx <= logLimX(lenAlpha + 1), SK_"@setPiwiPowetoCDF(): The condition `logLimX(1) <= logx and. logx < logLimX(size(logLimX))` must hold. logLimX(1), logx, logLimX(size(logLimX)) = "//getStr([logLimX(1), logx, logLimX(lenAlpha + 1)]))
170 9176 : CHECK_ASSERTION(__LINE__, isAscending(logLimX), SK_"@setPiwiPowetoCDF(): The conditions `isAscending(logLimX)` must hold. logLimX = "//getStr(logLimX))
171 : !CHECK_ASSERTION(__LINE__, logLimX(size(logLimX,1,IK)) < log(huge(logLimX)), SK_"@getPiwiPowetoLogPDFNF(): The condition `logLimX(size(logLimX)) < log(huge(logLimX))` must hold. logLimX = "//getStr(logLimX))
172 : !CHECK_ASSERTION(__LINE__, merge(alpha(lenAlpha) < -1._RKC, .true., size(logLimX,1,IK) == lenAlpha), SK_"@getPiwiPowetoLogPDFNF(): The condition `alpha(size(alpha)) < -1._RKC` must hold. alpha = "//getStr(alpha))
173 25136 : do i = lenAlpha, 1_IK, -1_IK
174 25136 : if (logx < logLimX(i)) cycle
175 9176 : call setPiwiPowetoCDF(cdf, logx, alpha, logLimX, logPDFNF, cumSumArea, i)
176 25136 : return
177 : end do
178 :
179 : #else
180 : !%%%%%%%%%%%%%%%%%%%%%%%%
181 : #error "Unrecognized interface."
182 : !%%%%%%%%%%%%%%%%%%%%%%%%
183 : #endif
|
