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_distGeomCyclic](@ref pm_distGeomCyclic).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Monday March 6, 2017, 3:22 pm, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin.<br>
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
28 : #if getGeomCyclicLogPMF_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 1074 : CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@getGeomCyclicLogPMF(): The condition `0 <= probSuccess` must hold. probSuccess = "//getStr(probSuccess))
32 1074 : CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@getGeomCyclicLogPMF(): The condition `probSuccess <= 1` must hold. probSuccess = "//getStr(probSuccess))
33 1074 : call setGeomCyclicLogPMF(logPMF, stepSuccess, log(probSuccess), period)
34 :
35 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
36 : #elif setGeomCyclicLogPMF_ENABLED
37 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
38 :
39 : real(RKC), parameter :: logProbSuccessMin = log(epsilon(0._RKC))
40 : real(RKC), parameter :: NEGINF = -log(huge(logPMF))
41 : real(RKC) :: logDenominator
42 : #if Def_ENABLED
43 : real(RKC) :: logProbFailure
44 : #endif
45 : ! Validate the input.
46 1146 : CHECK_ASSERTION(__LINE__, all(0_IK < [stepSuccess]), SK_"@setGeomCyclicLogPMF(): The condition `all(0 < [stepSuccess])` must hold. stepSuccess = "//getStr(stepSuccess))
47 3436 : CHECK_ASSERTION(__LINE__, all([stepSuccess] <= period), SK_"@setGeomCyclicLogPMF(): The condition `all([stepSuccess] <= period)` must hold. stepSuccess, period = "//getStr([real(RKC) :: stepSuccess, period]))
48 1123 : CHECK_ASSERTION(__LINE__, logProbSuccess <= 0._RKC, SK_"@setGeomCyclicLogPMF(): The condition `logProbSuccess <= 0.` must hold. logProbSuccess = "//getStr(logProbSuccess))
49 : ! Compute the PMF.
50 1123 : if (logProbSuccess < 0._RKC) then ! imperfect probability of success.
51 1121 : if (logProbSuccessMin < logProbSuccess) then
52 : #if Log_ENABLED
53 6 : CHECK_ASSERTION(__LINE__, abs(1._RKC - exp(logProbSuccess) - exp(logProbFailure)) < epsilon(0._RKC) * 100, SK_"@setGeomCyclicLogPMF(): The condition `exp(logProbFailure) + exp(logProbSuccess) == 1.` must hold. logProbFailure, logProbSuccess = "//getStr([logProbFailure, logProbSuccess]))
54 : #elif Def_ENABLED
55 1119 : logProbFailure = log(get1mexp(logProbSuccess))
56 : #else
57 : #error "Unrecognized interface."
58 : #endif
59 1121 : logDenominator = log(get1mexp(period * logProbFailure))
60 1144 : logPMF = logProbSuccess + (stepSuccess - 1_IK) * logProbFailure - logDenominator
61 : else
62 0 : logPMF = NEGINF
63 : end if
64 : else ! 100% probability of success.
65 : #if D0_ENABLED
66 2 : if (1_IK < stepSuccess) then ! stepSuccess can be larger than 1.
67 2 : logPMF = NEGINF
68 : else
69 0 : logPMF = 0._RKC
70 : end if
71 : #elif D1_ENABLED
72 : block
73 : integer(IK) :: istep
74 0 : do concurrent(istep = 1_IK : size(stepSuccess, 1, IK))
75 0 : if (1_IK < stepSuccess(istep)) then
76 0 : logPMF(istep) = -log(huge(logPMF))
77 : else
78 0 : logPMF(istep) = 0._RKC
79 : end if
80 : end do
81 : end block
82 : #endif
83 : end if
84 :
85 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
86 : #elif getGeomCyclicLogCDF_ENABLED
87 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
88 :
89 50 : CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@setGeomCyclicLogCDF(): The condition `0 <= probSuccess` must hold. probSuccess = "//getStr(probSuccess))
90 50 : CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@setGeomCyclicLogCDF(): The condition `probSuccess <= 1` must hold. probSuccess = "//getStr(probSuccess))
91 50 : call setGeomCyclicLogCDF(logCDF, stepSuccess, log(probSuccess), period)
92 :
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
94 : #elif setGeomCyclicLogCDF_ENABLED
95 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
96 :
97 : real(RKC) :: logDenominator
98 : #if Def_ENABLED
99 : real(RKC) :: logProbFailure
100 : #endif
101 : ! Validate the input.
102 101 : CHECK_ASSERTION(__LINE__, all(0_IK < [stepSuccess]), SK_"@setGeomCyclicLogPMF(): The condition `all(0 < [stepSuccess])` must hold. stepSuccess = "//getStr(stepSuccess))
103 302 : CHECK_ASSERTION(__LINE__, all([stepSuccess] <= period), SK_"@setGeomCyclicLogPMF(): The condition `all([stepSuccess] <= period)` must hold. stepSuccess, period = "//getStr([real(RKC) :: stepSuccess, period]))
104 98 : CHECK_ASSERTION(__LINE__, logProbSuccess <= 0._RKC, SK_"@setGeomCyclicLogPMF(): The condition `logProbSuccess <= 0.` must hold. logProbSuccess = "//getStr(logProbSuccess))
105 : ! Compute the CDF.
106 98 : if (logProbSuccess < 0._RKC) then ! imperfect probability of success.
107 : #if Log_ENABLED
108 3 : CHECK_ASSERTION(__LINE__, abs(1._RKC - exp(logProbSuccess) - exp(logProbFailure)) < epsilon(0._RKC) * 100, SK_"@setGeomCyclicLogPMF(): The condition `exp(logProbFailure) + exp(logProbSuccess) == 1.` must hold. logProbFailure, logProbSuccess = "//getStr([logProbFailure, logProbSuccess]))
109 : #elif Def_ENABLED
110 95 : logProbFailure = log(get1mexp(logProbSuccess))
111 : #else
112 : #error "Unrecognized interface."
113 : #endif
114 96 : logDenominator = log(get1mexp(logProbFailure * period))
115 102 : logCDF = log(get1mexp(logProbFailure * stepSuccess)) - logDenominator
116 : else ! 100% probability of success.
117 2 : logCDF = 0._RKC
118 : end if
119 :
120 : !%%%%%%%%%%%%%%%%%%%%%%%%
121 : #elif getGeomCyclicRand_ENABLED
122 : !%%%%%%%%%%%%%%%%%%%%%%%%
123 :
124 22006 : CHECK_ASSERTION(__LINE__, 0._RKC < probSuccess, SK_"@getGeomCyclicRand(): The condition `0 <= probSuccess` must hold. probSuccess = "//getStr(probSuccess))
125 22006 : if (probSuccess < 1._RKC) then
126 22005 : call setGeomCyclicRand(rand, log(1._RKC - probSuccess), period)
127 : else
128 1 : rand = 1_IK
129 1 : CHECK_ASSERTION(__LINE__, 0_IK < period, SK_"@getGeomCyclicRand(): The condition `0 < period` must hold. probSuccess = "//getStr(period))
130 1 : CHECK_ASSERTION(__LINE__, probSuccess <= 1._RKC, SK_"@getGeomCyclicRand(): The condition `probSuccess <= 1` must hold. probSuccess = "//getStr(probSuccess))
131 : end if
132 :
133 : !%%%%%%%%%%%%%%%%%%%%%%%%
134 : #elif setGeomCyclicRand_ENABLED
135 : !%%%%%%%%%%%%%%%%%%%%%%%%
136 :
137 42011 : if (1_IK < period) then
138 : #if RNGD_ENABLED
139 42010 : call setGeomRand(rand, logProbFailure)
140 : #elif RNGF_ENABLED || RNGX_ENABLED
141 1 : call setGeomRand(rng, rand, logProbFailure)
142 : #else
143 : #error "Unrecognized interface."
144 : #endif
145 : #if D0_ENABLED
146 42009 : rand = mod(rand, period)
147 42009 : if (rand == 0_IK) rand = period
148 : #elif D1_ENABLED
149 : block
150 : integer(IK) :: i
151 2 : do concurrent(i = 1 : size(rand, 1, IK))
152 4 : rand(i) = mod(rand(i), period)
153 6 : if (rand(i) == 0_IK) rand(i) = period
154 : end do
155 : end block
156 : #endif
157 : else
158 0 : rand = 1_IK
159 0 : CHECK_ASSERTION(__LINE__, 0_IK < period, SK_"@setGeomCyclicRand(): The condition `0 < period` must hold. probSuccess = "//getStr(period))
160 : end if
161 :
162 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163 : #elif isFailedGeomCyclicFit_ENABLED
164 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
165 :
166 : integer(IK) :: niter
167 : real(RKC) :: sumFreq
168 : real(RKC) :: logProbFailure
169 : real(RKC) :: logProbSuccess
170 2 : real(RKC) :: diff(size(stepSuccess, 1, IK))
171 1 : real(RKC) :: logFreqSuccess(size(freqSuccess, 1, IK))
172 : real(RKC), parameter :: small = sqrt(epsilon(0._RKC))
173 : integer(IK), parameter :: niter_def = int(1000 * precision(0._RKC) / 53., IK)
174 21 : CHECK_ASSERTION(__LINE__, all(0_IK < freqSuccess), SK_"@isFailedGeomCyclicFit(): The condition `all(0 < freqSuccess)` must hold. freqSuccess = "//getStr(freqSuccess))
175 21 : CHECK_ASSERTION(__LINE__, all(0_IK < stepSuccess), SK_"@isFailedGeomCyclicFit(): The condition `all(0 < stepSuccess)` must hold. stepSuccess = "//getStr(stepSuccess))
176 1 : CHECK_ASSERTION(__LINE__, 1_IK < size(stepSuccess, 1, IK), SK_"@isFailedGeomCyclicFit(): The condition `1 < size(stepSuccess)` must hold. size(stepSuccess) = "//getStr(size(stepSuccess, 1, IK)))
177 43 : CHECK_ASSERTION(__LINE__, maxval(stepSuccess) <= period, SK_"@isFailedGeomCyclicFit(): The condition `maxval(stepSuccess) <= period` must hold. maxval(stepSuccess), period = "//getStr([real(RKC) :: maxval(stepSuccess), period]))
178 3 : CHECK_ASSERTION(__LINE__, size(stepSuccess, 1, IK) <= period, SK_"@isFailedGeomCyclicFit(): The condition `size(stepSuccess) <= period` must hold. size(stepSuccess), period = "//getStr([size(stepSuccess, 1, IK), period]))
179 3 : CHECK_ASSERTION(__LINE__, size(stepSuccess, 1, IK) == size(freqSuccess, 1, IK), SK_"@isFailedGeomCyclicFit(): The condition `size(stepSuccess) == size(freqSuccess)` must hold. size(stepSuccess), size(freqSuccess) = "//getStr([size(stepSuccess, 1, IK), size(freqSuccess, 1, IK)]))
180 :
181 : ! First find the max likelihood.
182 :
183 1 : niter = niter_def
184 21 : sumFreq = sum(freqSuccess)
185 : !block
186 : ! real(RKC) :: xmin(1)
187 : ! xmin = .5_RKC
188 : ! failed = isFailedMinPowell(getNegLogLike, xmin)
189 : ! probSuccess = xmin(1)
190 : !end block
191 1 : probSuccess = getMinBrent(getNegLogLike, xlow = small, xupp = 1._RKC - small, niter = niter)
192 1 : failed = niter_def < niter
193 1 : if (failed) return
194 :
195 : ! Find the normalization factor.
196 :
197 1 : logProbFailure = log(1._RKC - probSuccess)
198 1 : logProbSuccess = log(probSuccess)
199 :
200 1 : niter = niter_def
201 1 : call setGeomCyclicLogPMF(diff, stepSuccess, logProbSuccess, logProbFailure, period)
202 21 : logFreqSuccess = log(real(freqSuccess, RKC))
203 21 : diff = logFreqSuccess - diff
204 1 : normFac = exp(getMinBrent(getSumDistSq, niter = niter))
205 1 : failed = niter_def < niter
206 :
207 : contains
208 :
209 21 : PURE function getNegLogLike(probSuccess) result(negLogLike)
210 : real(RKC), intent(in) :: probSuccess!(1)
211 : real(RKC) :: negLogLike
212 21 : if (0._RKC < probSuccess .and. probSuccess < 1._RKC) then
213 : !negLogLike = -sum(getGeomCyclicLogPMF(stepSuccess, probSuccess, period))
214 21 : negLogLike = (1 - probSuccess)**period
215 21 : if (negLogLike < 1) then
216 441 : negLogLike = sumFreq * (log(1 - negLogLike) - log(probSuccess)) + sum(freqSuccess * (1 - stepSuccess)) * log(1 - probSuccess)
217 : else
218 0 : negLogLike = sqrt(huge(0._RKC))
219 : end if
220 : else
221 0 : error stop "This is a miracle! Please inform the library developers at: https://github.com/cdslaborg/paramonte"
222 : !negLogLike = sqrt(huge(0._RKC))
223 : end if
224 21 : end function
225 :
226 12 : PURE function getSumDistSq(logNormFac) result(sumDistSq)
227 : real(RKC), intent(in) :: logNormFac
228 : real(RKC) :: sumDistSq
229 252 : sumDistSq = sum((diff - logNormFac)**2)
230 12 : end function
231 :
232 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
233 : #elif isFailedGeomCyclicFit_ENABLED && Def_ENABLED && 0
234 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
235 :
236 : integer(IK) :: i, counter
237 : real(RKC) :: logFreqSuccess(size(freqSuccess, 1, IK))
238 : integer(IK) :: stepSuccess(size(freqSuccess, 1, IK))
239 : CHECK_ASSERTION(__LINE__, size(stepSuccess, 1, IK) == size(freqSuccess, 1, IK), SK_"@isFailedGeomCyclicFit(): The condition `size(stepSuccess) == size(freqSuccess)` must hold. size(stepSuccess), size(freqSuccess) = "//getStr([size(stepSuccess, 1, IK), size(freqSuccess, 1, IK)]))
240 : counter = 0_IK
241 : do i = 1, size(stepSuccess, 1, IK)
242 : if (0 < freqSuccess(i)) then
243 : counter = counter + 1_IK
244 : logFreqSuccess(counter) = log(real(freqSuccess(counter), RKC))
245 : stepSuccess(counter) = stepSuccess(i)
246 : else
247 : CHECK_ASSERTION(__LINE__, 0 == freqSuccess(i), SK_"@isFailedGeomCyclicFit(): The condition `all(0 <= freqSuccess)` must hold. i, freqSuccess(i) = "//getStr([i, freqSuccess(i)]))
248 : end if
249 : end do
250 : failed = isFailedGeomCyclicFit(stepSuccess(1 : counter), logFreqSuccess(1 : counter), period, probSuccess, normFac)
251 :
252 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 : #elif isFailedGeomCyclicFit_ENABLED && Per_ENABLED && 0
254 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255 :
256 : real(RKC) :: successProbFisherTransLogNormFac(2) ! vector of the two parameters to fit.
257 : real(RKC), parameter :: SUCCESS_PROB_INIT_GUESS = 0.23_RKC
258 : real(RKC), parameter :: FISHER_TRANS_SUCCESS_PROB_INIT_GUESS = atanh(2 * (SUCCESS_PROB_INIT_GUESS - 0.5_RKC))
259 : CHECK_ASSERTION(__LINE__, all(0_IK < stepSuccess), SK_"@isFailedGeomCyclicFit(): The condition `all(0 < stepSuccess)` must hold. stepSuccess = "//getStr(stepSuccess))
260 : CHECK_ASSERTION(__LINE__, all(0._RKC <= logFreqSuccess), SK_"@isFailedGeomCyclicFit(): The condition `all(0 < stepSuccess)` must hold. stepSuccess = "//getStr(logFreqSuccess))
261 : CHECK_ASSERTION(__LINE__, 1_IK < size(stepSuccess, 1, IK), SK_"@isFailedGeomCyclicFit(): The condition `1 < size(stepSuccess)` must hold. size(stepSuccess) = "//getStr(size(stepSuccess, 1, IK)))
262 : CHECK_ASSERTION(__LINE__, maxval(stepSuccess) <= period, SK_"@isFailedGeomCyclicFit(): The condition `maxval(stepSuccess) <= period` must hold. maxval(stepSuccess), period = "//getStr([real(RKC) :: maxval(stepSuccess), period]))
263 : CHECK_ASSERTION(__LINE__, size(stepSuccess, 1, IK) <= period, SK_"@isFailedGeomCyclicFit(): The condition `size(stepSuccess) <= period` must hold. size(stepSuccess), period = "//getStr([size(stepSuccess, 1, IK), period]))
264 : CHECK_ASSERTION(__LINE__, size(stepSuccess, 1, IK) == size(logFreqSuccess, 1, IK), SK_"@isFailedGeomCyclicFit(): The condition `size(stepSuccess) == size(logFreqSuccess)` must hold. size(stepSuccess), size(logFreqSuccess) = "//getStr([size(stepSuccess, 1, IK), size(logFreqSuccess, 1, IK)]))
265 :
266 : ! Do Fisher transformation to ensure stability.
267 : successProbFisherTransLogNormFac(2) = 0._RKC
268 : successProbFisherTransLogNormFac(1) = FISHER_TRANS_SUCCESS_PROB_INIT_GUESS
269 : failed = isFailedMinPowell(getSumDistSq, successProbFisherTransLogNormFac)
270 : if (.not. failed) then
271 : probSuccess = 0.5_RKC * tanh(successProbFisherTransLogNormFac(1)) + 0.5_RKC ! reverse Fisher-transform
272 : normFac = exp(successProbFisherTransLogNormFac(2))
273 : end if
274 :
275 : contains
276 :
277 : PURE function getSumDistSq(successProbFisherTransLogNormFac) result(sumDistSq)
278 : real(RKC), intent(in) :: successProbFisherTransLogNormFac(2)
279 : real(RKC) :: logPMF(size(stepSuccess, 1, IK)), logProbSuccess, sumDistSq
280 : logProbSuccess = log(0.5_RKC * tanh(successProbFisherTransLogNormFac(1)) + 0.5_RKC)
281 : call setGeomCyclicLogPMF(logPMF, stepSuccess, logProbSuccess, period)
282 : sumDistSq = sum((logFreqSuccess - logPMF - successProbFisherTransLogNormFac(2))**2)
283 : end function getSumDistSq
284 : #else
285 : !%%%%%%%%%%%%%%%%%%%%%%%%
286 : #error "Unrecognized interface."
287 : !%%%%%%%%%%%%%%%%%%%%%%%%
288 : #endif
|
