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 the tests of [pm_distPower](@ref pm_distPower).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, 12:27 AM Tuesday, February 22, 2022, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if CK_ENABLED
28 : use pm_complexAbs, only: abs, log, operator(<), operator(<=)
29 : #elif !RK_ENABLED
30 : #error "Unrecognized interface."
31 : #endif
32 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 : #if getPowerLogPDF_ENABLED || setPowerLogPDF_ENABLED
34 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 :
36 : integer(IK) :: i
37 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
38 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
39 : real(RKC) , parameter :: TOL = sqrt(epsilon(0._RKC))
40 : real(RKC) :: alpha, logPDFNF, logMinX, logMaxX
41 : real(RKC) :: output, output_ref, diff
42 :
43 8 : assertion = .true._LK
44 :
45 408 : do i = 1_IK, 50_IK
46 400 : logMinX = getUnifRand(-5._RKC, +4._RKC)
47 400 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
48 400 : call runTestsWith(logMinX)
49 408 : call runTestsWith()
50 : end do
51 :
52 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53 :
54 : contains
55 :
56 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 :
58 800 : subroutine runTestsWith(logMinX)
59 :
60 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(log(huge(0._RKC)))
61 : real(RKC), intent(in), optional :: logMinX
62 :
63 800 : alpha = getUnifRand(0.5_RKC, +3._RKC)
64 800 : if (present(logMinX)) then
65 400 : logPDFNF = getPowerLogPDFNF(alpha, logMinX, logMaxX)
66 : else
67 400 : logPDFNF = getPowerLogPDFNF(alpha, logMaxX)
68 : end if
69 :
70 : #if getPowerLogPDF_ENABLED
71 : block
72 : real(RKC) :: logx
73 400 : logx = getUnifRand(getOption(-SQRT_LOG_HUGE, logMinX), logMaxX)
74 400 : call setPowerLogPDF(output_ref, logx, alpha, logPDFNF)
75 400 : if (present(logMinX)) then
76 200 : output = getPowerLogPDF(logx, alpha, logMinX, logMaxX)
77 : else
78 200 : output = getPowerLogPDF(logx, alpha, logMaxX)
79 : end if
80 400 : call compare(logMinX, logx)
81 : end block
82 : #elif setPowerLogPDF_ENABLED
83 400 : output_ref = 1._RKC
84 : block
85 : use pm_quadPack, only: isFailedQuad, getQuadErr, GK31, weps
86 : integer(IK) :: err, neval, nint, sindex(2000)
87 : real(RKC) :: lb, abserr, sinfo(4,2000)
88 : character(255, SK) :: msg
89 400 : msg = SK_" "
90 400 : lb = 0._RKC
91 400 : if (present(logMinX)) lb = exp(logMinX)
92 400 : if (alpha > 1._RKC) then
93 318 : assertion = assertion .and. .not. isFailedQuad(getPowerPDF, lb, exp(logMaxX), output, reltol = TOL, msg = msg)
94 : else
95 82 : err = getQuadErr(getPowerPDF, lb, exp(logMaxX), 0._RKC, TOL, GK31, weps, output, abserr, sinfo, sindex, neval, nint) ! extends QAGS/QAGI routines of QuadPack.
96 82 : assertion = assertion .and. err == 0_IK
97 : end if
98 400 : call report(logMinX)
99 400 : call test%assert(assertion, SK_"@setPowerLogPDF(): The integral of the Power distribution must be computed without failure. msg: "//trim(msg), int(__LINE__, IK))
100 : end block
101 : call compare()
102 400 : call report(logMinX)
103 400 : call test%assert(assertion, SK_"@setPowerLogPDF(): The PDF of the Power distribution must be computed correctly.", int(__LINE__, IK))
104 : #else
105 : #error "Unrecognized interface."
106 : #endif
107 800 : end subroutine
108 :
109 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110 :
111 : #if getPowerLogPDF_ENABLED
112 400 : subroutine compare(logMinX, logx)
113 : real(RKC), intent(in), optional :: logMinX, logx
114 400 : diff = abs(output - output_ref)
115 400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
116 400 : call report(logMinX, logx)
117 400 : call test%assert(assertion, SK_"@getPowerLogPDF(): The PDF of the Power distribution must be computed correctly.", int(__LINE__, IK))
118 400 : end subroutine
119 : #elif setPowerLogPDF_ENABLED
120 : subroutine compare()
121 400 : diff = abs(output - output_ref)
122 400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
123 : end subroutine
124 106934 : function getPowerPDF(x) result(pdf)
125 : real(RKC), intent(in) :: x
126 : real(RKC) :: pdf
127 106934 : call setPowerLogPDF(pdf, log(x), alpha, logPDFNF)
128 106934 : pdf = exp(pdf)
129 106934 : end function
130 : #else
131 : #error "Unrecognized interface."
132 : #endif
133 :
134 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135 :
136 1200 : subroutine report(logMinX, logx)
137 : real(RKC), intent(in), optional :: logMinX, logx
138 1200 : if (test%traceable .and. .not. assertion) then
139 : ! LCOV_EXCL_START
140 : write(test%disp%unit,"(*(g0,:,', '))")
141 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
142 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMinX)", present(logMinX)
143 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
144 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMinX)", getOption(-HUGE_RKC, logMinX)
145 : write(test%disp%unit,"(*(g0,:,', '))") "logPDFNF", logPDFNF
146 : write(test%disp%unit,"(*(g0,:,', '))") "logMaxX", logMaxX
147 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
148 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
149 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
150 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
151 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
152 : write(test%disp%unit,"(*(g0,:,', '))")
153 : ! LCOV_EXCL_STOP
154 : end if
155 1200 : end subroutine
156 :
157 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158 :
159 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160 : #elif getPowerLogCDF_ENABLED || setPowerLogCDF_ENABLED
161 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162 :
163 : integer(IK) :: i
164 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
165 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
166 : real(RKC) , parameter :: TOL = sqrt(epsilon(0._RKC))
167 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
168 : real(RKC) :: output, output_ref, diff
169 : logical(LK) :: logMinXPresent
170 :
171 8 : assertion = .true._LK
172 :
173 2408 : do i = 1_IK, 300_IK
174 2400 : logMinX = getUnifRand(-5._RKC, +4._RKC)
175 2400 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
176 2400 : call runTestsWith(logMinX = logMinX)
177 2408 : call runTestsWith()
178 : end do
179 :
180 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
181 :
182 : contains
183 :
184 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185 :
186 4800 : subroutine runTestsWith(logMinX)
187 :
188 : use pm_val2str, only: getStr
189 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
190 : real(RKC), intent(in), optional :: logMinX
191 : real(RKC) :: logx
192 :
193 4800 : logMinXPresent = present(logMinX)
194 4800 : alpha = getUnifRand(0.5_RKC, +3._RKC)
195 4800 : if (logMinXPresent) then
196 2400 : logCDFNF = getPowerLogCDFNF(alpha, logMinX, logMaxX)
197 : else
198 2400 : logCDFNF = getPowerLogCDFNF(alpha, logMaxX)
199 : end if
200 :
201 : #if getPowerLogCDF_ENABLED
202 :
203 2400 : logx = getUnifRand(getOption(-SQRT_LOG_HUGE, logMinX), logMaxX)
204 :
205 2400 : if (logMinXPresent) then
206 1200 : call setPowerLogCDF(output_ref, logx, alpha, logMinX, logCDFNF)
207 1200 : output = getPowerLogCDF(logx, alpha, logMinX, logMaxX)
208 : else
209 1200 : call setPowerLogCDF(output_ref, logx, alpha, logCDFNF)
210 1200 : output = getPowerLogCDF(logx, alpha, logMaxX)
211 : end if
212 2400 : call compare(__LINE__, logMinX, logx)
213 :
214 : #elif setPowerLogCDF_ENABLED
215 :
216 2400 : output_ref = -LOG_HUGE
217 2400 : logx = getOption(-LOG_HUGE, logMinX)
218 :
219 2400 : if (logMinXPresent) then
220 1200 : call setPowerLogCDF(output, logx, alpha, logMinX, logCDFNF)
221 : else
222 1200 : call setPowerLogCDF(output, logx, alpha, logCDFNF)
223 : end if
224 2400 : call report(logMinX, logx)
225 2400 : call test%assert(assertion, SK_"@setPowerLogCDF(): The CDF of the Power distribution at the lower limit of the support must be zero.", int(__LINE__, IK))
226 :
227 2400 : output_ref = 0._RKC
228 2400 : logx = logMaxX
229 :
230 2400 : if (logMinXPresent) then
231 1200 : call setPowerLogCDF(output, logx, alpha, logMinX, logCDFNF)
232 : else
233 1200 : call setPowerLogCDF(output, logx, alpha, logCDFNF)
234 : end if
235 2400 : call report(logMinX, logx)
236 2400 : call test%assert(assertion, SK_"@setPowerLogCDF(): The CDF of the Power distribution at the upper limit of the support must be unity.", int(__LINE__, IK))
237 :
238 : ! Compare the CDF difference at random points with the corresponding integral of the PDF.
239 :
240 : block
241 : use pm_quadPack, only: isFailedQuad, getQuadErr, GK31, weps
242 : integer(IK) :: err, neval, nint, sindex(2000)
243 : real(RKC) :: lb, abserr, sinfo(4,2000)
244 : character(255, SK) :: msg
245 2400 : msg = SK_" "
246 2400 : lb = 0._RKC
247 2400 : if (logMinXPresent) lb = exp(logMinX)
248 2400 : logx = log(getUnifRand((exp(logMaxX) + lb) * .5_RKC, exp(logMaxX)))
249 2400 : if (alpha > 1._RKC) then
250 1925 : assertion = assertion .and. .not. isFailedQuad(getPowerPDF, lb, exp(logx), output_ref, reltol = TOL, msg = msg)
251 : else
252 475 : err = getQuadErr(getPowerPDF, lb, exp(logx), 0._RKC, TOL, GK31, weps, output_ref, abserr, sinfo, sindex, neval, nint) ! extends QAGS/QAGI routines of QuadPack.
253 475 : assertion = assertion .and. err == 0_IK
254 : end if
255 2400 : call report(logMinX, logx)
256 2400 : call test%assert(assertion, SK_"@setPowerLogCDF(): The integral of the Power distribution must be computed without failure. msg: "//trim(msg), int(__LINE__, IK))
257 2400 : if (logMinXPresent) then
258 1200 : call setPowerLogCDF(output, logx, alpha, logMinX, logCDFNF)
259 : else
260 1200 : call setPowerLogCDF(output, logx, alpha, logCDFNF)
261 : end if
262 2400 : output = exp(output)
263 : call compare()
264 2400 : call report(logMinX, logx)
265 2400 : call test%assert(assertion, SK_"@setPowerLogCDF(): The CDF of the Power distribution must be computed correctly.", int(__LINE__, IK))
266 : end block
267 :
268 : #else
269 : #error "Unrecognized interface."
270 : #endif
271 4800 : end subroutine
272 :
273 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
274 :
275 : #if getPowerLogCDF_ENABLED
276 2400 : subroutine compare(line, logMinX, logx)
277 : integer :: line
278 : real(RKC), intent(in), optional :: logMinX, logx
279 2400 : diff = abs(output - output_ref)
280 2400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
281 2400 : call report(logMinX, logx)
282 2400 : call test%assert(assertion, SK_"@getPowerLogCDF(): The CDF of the Power distribution must be computed correctly.", int(line, IK))
283 2400 : end subroutine
284 : #elif setPowerLogCDF_ENABLED
285 : subroutine compare()
286 2400 : diff = abs(output - output_ref)
287 2400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
288 : end subroutine
289 639862 : function getPowerPDF(x) result(pdf)
290 : use pm_distPower, only: getPowerLogPDF
291 : real(RKC), intent(in) :: x
292 : real(RKC) :: pdf
293 639862 : if (logMinXPresent) then
294 65514 : pdf = exp(getPowerLogPDF(log(x), alpha, logMinX, logMaxX))
295 : else
296 574348 : pdf = exp(getPowerLogPDF(log(x), alpha, logMaxX))
297 : end if
298 639862 : end function
299 : #else
300 : #error "Unrecognized interface."
301 : #endif
302 :
303 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304 :
305 12000 : subroutine report(logMinX, logx)
306 : real(RKC), intent(in), optional :: logMinX, logx
307 12000 : if (test%traceable .and. .not. assertion) then
308 : ! LCOV_EXCL_START
309 : write(test%disp%unit,"(*(g0,:,', '))")
310 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
311 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMinX)", present(logMinX)
312 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
313 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMinX)", getOption(-HUGE_RKC, logMinX)
314 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
315 : write(test%disp%unit,"(*(g0,:,', '))") "logMaxX", logMaxX
316 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
317 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
318 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
319 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
320 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
321 : write(test%disp%unit,"(*(g0,:,', '))")
322 : ! LCOV_EXCL_STOP
323 : end if
324 12000 : end subroutine
325 :
326 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 :
328 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
329 : #elif getPowerLogQuan_ENABLED || setPowerLogQuan_ENABLED
330 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
331 :
332 : integer(IK) :: i
333 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
334 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
335 : real(RKC) , parameter :: TOL = sqrt(epsilon(0._RKC))! * 1000
336 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
337 : real(RKC) :: logx, output, output_ref, diff
338 : logical(LK) :: logMinXPresent
339 :
340 8 : assertion = .true._LK
341 :
342 4008 : do i = 1_IK, 500_IK
343 4000 : logMinX = getUnifRand(-LOG_HUGE, LOG_HUGE)
344 4000 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
345 4000 : call runTestsWith(logMinX = logMinX)
346 4008 : call runTestsWith()
347 : end do
348 :
349 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
350 :
351 : contains
352 :
353 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 :
355 8000 : subroutine runTestsWith(logMinX)
356 :
357 : use pm_val2str, only: getStr
358 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
359 : real(RKC), intent(in), optional :: logMinX
360 :
361 8000 : logMinXPresent = present(logMinX)
362 8000 : alpha = getUnifRand(0.5_RKC, +3._RKC)
363 8000 : if (logMinXPresent) then
364 4000 : logCDFNF = getPowerLogCDFNF(alpha, logMinX, logMaxX)
365 : else
366 4000 : logCDFNF = getPowerLogCDFNF(alpha, logMaxX)
367 : end if
368 :
369 8000 : if (i > 1_IK) then
370 7984 : logx = getUnifRand(getOption(logMaxX - 4._RKC, logMinX), logMaxX)
371 : else ! test for logCDF = 0.
372 16 : logx = logMaxX
373 16 : output_ref = 0._RKC ! logCDF
374 : end if
375 8000 : if (logMinXPresent) then
376 4000 : if (i > 1_IK) call setPowerLogCDF(output_ref, logx, alpha, logMinX, logCDFNF)
377 : #if getPowerLogQuan_ENABLED
378 2000 : output = getPowerLogQuan(output_ref, alpha, logMinX, logMaxX)
379 : #elif setPowerLogQuan_ENABLED
380 2000 : call setPowerLogQuan(output, output_ref, alpha, logMinX, logCDFNF)
381 : #else
382 : #error "Unrecognized interface."
383 : #endif
384 : else
385 4000 : if (i > 1_IK) call setPowerLogCDF(output_ref, logx, alpha, logCDFNF)
386 : #if getPowerLogQuan_ENABLED
387 2000 : output = getPowerLogQuan(output_ref, alpha, logMaxX)
388 : #elif setPowerLogQuan_ENABLED
389 2000 : call setPowerLogQuan(output, output_ref, alpha, logCDFNF)
390 : #endif
391 : end if
392 8000 : output_ref = logx
393 8000 : call report(__LINE__, logMinX, logx)
394 :
395 8000 : end subroutine
396 :
397 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
398 :
399 8000 : subroutine report(line, logMinX, logx)
400 : integer, intent(in) :: line
401 : real(RKC), intent(in), optional :: logMinX, logx
402 8000 : diff = abs(output - output_ref)
403 8000 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
404 8000 : if (test%traceable .and. .not. assertion) then
405 : ! LCOV_EXCL_START
406 : write(test%disp%unit,"(*(g0,:,', '))")
407 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
408 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMinX)", present(logMinX)
409 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
410 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMinX)", getOption(-HUGE_RKC, logMinX)
411 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
412 : write(test%disp%unit,"(*(g0,:,', '))") "logMaxX", logMaxX
413 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
414 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
415 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
416 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
417 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
418 : write(test%disp%unit,"(*(g0,:,', '))")
419 : ! LCOV_EXCL_STOP
420 : end if
421 8000 : call test%assert(assertion, SK_"@getPowerLogQuan(): The Quantile of the Power distribution must be computed correctly.", int(line, IK))
422 8000 : end subroutine
423 :
424 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425 :
426 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427 : #elif getPowerLogRand_ENABLED || setPowerLogRand_ENABLED
428 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
429 :
430 : integer(IK) :: i
431 :
432 : integer(IK) , parameter :: NSIM = 1000_IK
433 : real(RKC) , parameter :: MEAN_REF = 0.5_RKC
434 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
435 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
436 : real(RKC) , parameter :: TOL = 0.1_RKC
437 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
438 : real(RKC) :: LogRand(NSIM), CDF(NSIM), mean, diff
439 : logical(LK) :: logMinXPresent
440 :
441 8 : assertion = .true._LK
442 :
443 24 : do i = 1_IK, 2_IK
444 16 : logMinX = getUnifRand(-LOG_HUGE, LOG_HUGE)
445 16 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
446 16 : call runTestsWith(logMinX = logMinX)
447 24 : call runTestsWith()
448 : end do
449 :
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 :
452 : contains
453 :
454 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
455 :
456 32 : subroutine runTestsWith(logMinX)
457 :
458 : use pm_val2str, only: getStr
459 : use pm_sampleMean, only: getMean
460 : use pm_distNegExp, only: getNegExpRand
461 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
462 : real(RKC), intent(in), optional :: logMinX
463 : integer(IK) :: j
464 :
465 32 : logMinXPresent = present(logMinX)
466 32 : alpha = getUnifRand(0.5_RKC, +3._RKC)
467 32 : if (logMinXPresent) then
468 16 : logCDFNF = getPowerLogCDFNF(alpha, logMinX, logMaxX)
469 : else
470 16 : logCDFNF = getPowerLogCDFNF(alpha, logMaxX)
471 : end if
472 :
473 32032 : LogRand = getUnifRand(getOption(logMaxX - 4._RKC, logMinX), logMaxX, s1 = size(LogRand, kind = IK))
474 32 : if (logMinXPresent) then
475 16016 : do j = 1_IK, size(LogRand, kind = IK)
476 : #if getPowerLogRand_ENABLED
477 8000 : LogRand(j) = getPowerLogRand(alpha, logMinX, logMaxX) ! Truncated Power distribution.
478 : #elif setPowerLogRand_ENABLED
479 8000 : call setPowerLogRand(LogRand(j), getNegExpRand(1._RKC), alpha, logMinX, logCDFNF)
480 : #else
481 : #error "Unrecognized interface."
482 : #endif
483 16000 : call setPowerLogCDF(CDF(j), LogRand(j), alpha, logMinX, logCDFNF)
484 16016 : call report(__LINE__, logMinX, LogRand(j))
485 : end do
486 : else
487 16016 : do j = 1_IK, size(LogRand, kind = IK)
488 : #if getPowerLogRand_ENABLED
489 8000 : LogRand(j) = getPowerLogRand(alpha, logMaxX) ! Truncated Power distribution.
490 : #elif setPowerLogRand_ENABLED
491 8000 : call setPowerLogRand(LogRand(j), getNegExpRand(1._RKC), alpha, logCDFNF)
492 : #else
493 : #error "Unrecognized interface."
494 : #endif
495 16000 : call setPowerLogCDF(CDF(j), LogRand(j), alpha, logCDFNF)
496 16016 : call report(__LINE__, logMinX, LogRand(j))
497 : end do
498 : end if
499 32032 : mean = getMean(exp(CDF))
500 32 : call report(__LINE__, logMinX, mean = mean)
501 :
502 32 : end subroutine
503 :
504 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
505 :
506 32032 : subroutine report(line, logMinX, logRand, mean)
507 : integer, intent(in) :: line
508 : real(RKC), intent(in), optional :: logMinX, logRand, mean
509 32032 : if (present(mean) .and. present(logRand)) then
510 : error stop "Internal error occurred. Both `logRand` and `mean` input arguments cannot be present simultaneously." ! LCOV_EXCL_LINE
511 32032 : elseif (present(mean)) then
512 32 : diff = abs(mean - MEAN_REF)
513 32 : assertion = assertion .and. diff <= max(abs(MEAN_REF) * TOL, TOL)
514 32000 : elseif (present(logRand)) then
515 32000 : assertion = assertion .and. logical(logMaxX >= logRand, LK) .and. logical(logRand >= getOption(logMinX, -huge(logMinX)), LK)
516 : else
517 : error stop "Internal error occurred. Either `logRand` and `mean` input argument must be present." ! LCOV_EXCL_LINE
518 : end if
519 32032 : if (test%traceable .and. .not. assertion) then
520 : ! LCOV_EXCL_START
521 : write(test%disp%unit,"(*(g0,:,', '))")
522 : write(test%disp%unit,"(*(g0,:,', '))") "present(LogRand)", present(LogRand)
523 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMinX)", present(logMinX)
524 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, LogRand)", getOption(HUGE_RKC, LogRand)
525 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMinX)", getOption(-HUGE_RKC, logMinX)
526 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
527 : write(test%disp%unit,"(*(g0,:,', '))") "logMaxX", logMaxX
528 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
529 : if (present(mean)) then
530 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
531 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
532 : write(test%disp%unit,"(*(g0,:,', '))") "mean", mean
533 : write(test%disp%unit,"(*(g0,:,', '))") "MEAN_REF", MEAN_REF
534 : write(test%disp%unit,"(*(g0,:,', '))")
535 : end if
536 : ! LCOV_EXCL_STOP
537 : end if
538 32032 : if (present(mean)) then
539 32 : call test%assert(assertion, SK_"The generated random value must be Power-distributed. Due to the statistical nature of the computations, this test can occasional fail, in which case, rerunning the test may fix the failure.", int(line, IK))
540 32000 : elseif (present(logRand)) then
541 32000 : call test%assert(assertion, SK_"The generated Power random value must be within the support of the distribution.", int(line, IK))
542 : end if
543 32032 : end subroutine
544 :
545 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
546 :
547 : #else
548 : !%%%%%%%%%%%%%%%%%%%%%%%%
549 : #error "Unrecognized interface."
550 : !%%%%%%%%%%%%%%%%%%%%%%%%
551 : #endif
|
