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_distPareto](@ref pm_distPareto).
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 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 : #if getParetoLogPDF_ENABLED || setParetoLogPDF_ENABLED
35 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 :
37 : integer(IK) :: i
38 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
39 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
40 : real(RKC) , parameter :: TOL = sqrt(epsilon(0._RKC))
41 : real(RKC) :: logx, alpha, logPDFNF, logMinX, logMaxX
42 : real(RKC) :: output, output_ref, diff
43 :
44 8 : assertion = .true._LK
45 :
46 168 : do i = 1_IK, 20_IK
47 160 : logMinX = getUnifRand(-5._RKC, +4._RKC)
48 160 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
49 160 : call runTestsWith(logMaxX)
50 168 : call runTestsWith()
51 : end do
52 :
53 : contains
54 :
55 320 : subroutine runTestsWith(logMaxX)
56 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(log(huge(0._RKC)))
57 : real(RKC), intent(in), optional :: logMaxX
58 320 : alpha = -getUnifRand(0.5_RKC, +3._RKC)
59 320 : if (present(logMaxX)) then
60 160 : logPDFNF = getParetoLogPDFNF(alpha, logMinX, logMaxX)
61 : else
62 160 : logPDFNF = getParetoLogPDFNF(alpha, logMinX)
63 : end if
64 : #if getParetoLogPDF_ENABLED
65 160 : logx = getUnifRand(logMinX, getOption(SQRT_LOG_HUGE, logMaxX))
66 160 : call setParetoLogPDF(output_ref, logx, alpha, logPDFNF)
67 160 : if (present(logMaxX)) then
68 80 : output = getParetoLogPDF(logx, alpha, logMinX, logMaxX)
69 : else
70 80 : output = getParetoLogPDF(logx, alpha, logMinX)
71 : end if
72 160 : call compare(__LINE__, logMaxX, logx)
73 : #elif setParetoLogPDF_ENABLED
74 160 : output_ref = 1._RKC
75 : block
76 : character(255, SK) :: msg
77 160 : msg = SK_" "
78 160 : logx = huge(0._RKC)
79 160 : if (present(logMaxX)) logx = exp(logMaxX)
80 160 : assertion = assertion .and. .not. isFailedQuad(getParetoPDF, exp(logMinX), logx, output, reltol = TOL, msg = msg)
81 160 : if (.not. present(logMaxX)) logx = LOG_HUGE
82 160 : call report(logMaxX, logx)
83 160 : call test%assert(assertion, SK_"@setParetoLogPDF(): The integral of the Pareto distribution must be computed without failure. msg: "//trim(msg), int(__LINE__, IK))
84 : end block
85 : call compare()
86 160 : call report(logMaxX)
87 160 : call test%assert(assertion, SK_"@setParetoLogPDF(): The PDF of the Pareto distribution must be computed correctly.", int(__LINE__, IK))
88 : #else
89 : #error "Unrecognized interface."
90 : #endif
91 320 : end subroutine
92 :
93 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 :
95 : #if getParetoLogPDF_ENABLED
96 160 : subroutine compare(line, logMaxX, logx)
97 : integer, intent(in) :: line
98 : real(RKC), intent(in), optional :: logMaxX, logx
99 160 : diff = abs(output - output_ref)
100 160 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
101 160 : call report(logMaxX, logx)
102 160 : call test%assert(assertion, SK_"@getParetoLogPDF(): The PDF of the Pareto distribution must be computed correctly.", int(line, IK))
103 160 : end subroutine
104 : #elif setParetoLogPDF_ENABLED
105 : subroutine compare()
106 160 : diff = abs(output - output_ref)
107 160 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
108 : end subroutine
109 73878 : function getParetoPDF(x) result(pdf)
110 : real(RKC), intent(in) :: x
111 : real(RKC) :: pdf
112 73878 : call setParetoLogPDF(pdf, log(x), alpha, logPDFNF)
113 73878 : pdf = exp(pdf)
114 73878 : end function
115 : #else
116 : #error "Unrecognized interface."
117 : #endif
118 :
119 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120 :
121 480 : subroutine report(logMaxX, logx)
122 : real(RKC), intent(in), optional :: logMaxX, logx
123 480 : if (test%traceable .and. .not. assertion) then
124 : ! LCOV_EXCL_START
125 : write(test%disp%unit,"(*(g0,:,', '))")
126 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
127 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMaxX)", present(logMaxX)
128 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
129 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMaxX)", getOption(-HUGE_RKC, logMaxX)
130 : write(test%disp%unit,"(*(g0,:,', '))") "logPDFNF", logPDFNF
131 : write(test%disp%unit,"(*(g0,:,', '))") "logMinX", logMinX
132 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
133 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
134 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
135 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
136 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
137 : write(test%disp%unit,"(*(g0,:,', '))")
138 : ! LCOV_EXCL_STOP
139 : end if
140 480 : end subroutine
141 :
142 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143 :
144 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 : #elif getParetoLogCDF_ENABLED || setParetoLogCDF_ENABLED
146 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 :
148 : integer(IK) :: i
149 :
150 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
151 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
152 : real(RKC) , parameter :: TOL = epsilon(0._RKC) * 10000
153 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
154 : real(RKC) :: output, output_ref, diff
155 : logical(LK) :: logMaxXPresent
156 :
157 8 : assertion = .true._LK
158 :
159 2408 : do i = 1_IK, 300_IK
160 2400 : logMinX = getUnifRand(-5._RKC, +4._RKC)
161 2400 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
162 2400 : call runTestsWith(logMaxX = logMaxX)
163 2408 : call runTestsWith()
164 : end do
165 :
166 : contains
167 :
168 4800 : subroutine runTestsWith(logMaxX)
169 :
170 : use pm_val2str, only: getStr
171 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
172 : real(RKC), intent(in), optional :: logMaxX
173 : real(RKC) :: logx
174 4800 : logMaxXPresent = present(logMaxX)
175 4800 : alpha = -getUnifRand(sqrt(epsilon(0._RKC)), +3._RKC)
176 4800 : if (logMaxXPresent) then
177 2400 : logCDFNF = getParetoLogCDFNF(alpha, logMinX, logMaxX)
178 : else
179 2400 : logCDFNF = 0._RKC
180 : end if
181 :
182 : #if getParetoLogCDF_ENABLED
183 :
184 2400 : logx = getUnifRand(logMinX, getOption(logMinX + 4._RKC, logMaxX))
185 2400 : if (logMaxXPresent) then
186 1200 : call setParetoLogCDF(output_ref, logx, alpha, logMinX, logCDFNF)
187 1200 : output = getParetoLogCDF(logx, alpha, logMinX, logMaxX)
188 : else
189 1200 : call setParetoLogCDF(output_ref, logx, alpha, logMinX)
190 1200 : output = getParetoLogCDF(logx, alpha, logMinX)
191 : end if
192 2400 : call compare(__LINE__, logMaxX, logx)
193 :
194 : #elif setParetoLogCDF_ENABLED
195 :
196 2400 : output_ref = -LOG_HUGE
197 2400 : logx = logMinX
198 2400 : if (logMaxXPresent) then
199 1200 : call setParetoLogCDF(output, logx, alpha, logMinX, logCDFNF)
200 : else
201 1200 : call setParetoLogCDF(output, logx, alpha, logMinX)
202 : end if
203 2400 : call report(logMaxX, logx)
204 2400 : call test%assert(assertion, SK_"@setParetoLogCDF(): The CDF of the Pareto distribution at the lower limit of the support must be zero.", int(__LINE__, IK))
205 :
206 2400 : output_ref = 0._RKC
207 2400 : logx = getOption(+LOG_HUGE, logMaxX)
208 2400 : if (logMaxXPresent) then
209 1200 : call setParetoLogCDF(output, logx, alpha, logMinX, logCDFNF)
210 : else
211 1200 : call setParetoLogCDF(output, logx, alpha, logMinX)
212 : end if
213 2400 : call report(logMaxX, logx)
214 2400 : call test%assert(assertion, SK_"@setParetoLogCDF(): The CDF of the Pareto distribution at the upper limit of the support must be unity.", int(__LINE__, IK))
215 :
216 : ! Compare the CDF difference at random points with the corresponding integral of the PDF.
217 :
218 : block
219 : character(255, SK) :: msg
220 2400 : msg = SK_" "
221 2400 : logx = getUnifRand(logMinX, getOption(logMinX + 4._RKC, logMaxX))
222 2400 : assertion = assertion .and. .not. isFailedQuad(getParetoPDF, exp(logMinX), exp(logx), output_ref, reltol = TOL, msg = msg)
223 2400 : call report(logMaxX, logx)
224 2400 : call test%assert(assertion, SK_"@setParetoLogCDF(): The integral of the Pareto distribution must be computed without failure. msg: "//trim(msg), int(__LINE__, IK))
225 2400 : if (logMaxXPresent) then
226 1200 : call setParetoLogCDF(output, logx, alpha, logMinX, logCDFNF)
227 : else
228 1200 : call setParetoLogCDF(output, logx, alpha, logMinX)
229 : end if
230 2400 : output = exp(output)
231 : call compare()
232 2400 : call report(logMaxX, logx)
233 2400 : call test%assert(assertion, SK_"@setParetoLogCDF(): The CDF of the Pareto distribution must be computed correctly.", int(__LINE__, IK))
234 : end block
235 :
236 : #else
237 : #error "Unrecognized interface."
238 : #endif
239 4800 : end subroutine
240 :
241 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
242 :
243 : #if getParetoLogCDF_ENABLED
244 2400 : subroutine compare(line, logMaxX, logx)
245 : integer :: line
246 : real(RKC), intent(in), optional :: logMaxX, logx
247 2400 : diff = abs(output - output_ref)
248 2400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL) * 10
249 2400 : call report(logMaxX, logx)
250 2400 : call test%assert(assertion, SK_"@getParetoLogCDF(): The CDF of the Pareto distribution must be computed correctly.", int(line, IK))
251 2400 : end subroutine
252 : #elif setParetoLogCDF_ENABLED
253 : subroutine compare()
254 2400 : diff = abs(output - output_ref)
255 2400 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL) * 10
256 : end subroutine
257 363888 : function getParetoPDF(x) result(pdf)
258 : use pm_distPareto, only: getParetoLogPDF
259 : real(RKC), intent(in) :: x
260 : real(RKC) :: pdf
261 363888 : if (logMaxXPresent) then
262 141834 : pdf = exp(getParetoLogPDF(log(x), alpha, logMinX, logMaxX))
263 : else
264 222054 : pdf = exp(getParetoLogPDF(log(x), alpha, logMinX))
265 : end if
266 363888 : end function
267 : #else
268 : #error "Unrecognized interface."
269 : #endif
270 12000 : subroutine report(logMaxX, logx)
271 : real(RKC), intent(in), optional :: logMaxX, logx
272 12000 : if (test%traceable .and. .not. assertion) then
273 : ! LCOV_EXCL_START
274 : write(test%disp%unit,"(*(g0,:,', '))")
275 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
276 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMaxX)", present(logMaxX)
277 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
278 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMaxX)", getOption(-HUGE_RKC, logMaxX)
279 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
280 : write(test%disp%unit,"(*(g0,:,', '))") "logMinX", logMinX
281 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
282 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
283 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
284 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
285 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
286 : write(test%disp%unit,"(*(g0,:,', '))")
287 : ! LCOV_EXCL_STOP
288 : end if
289 12000 : end subroutine
290 :
291 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292 : #elif getParetoLogQuan_ENABLED || setParetoLogQuan_ENABLED
293 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294 :
295 : integer(IK) :: i
296 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
297 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
298 : real(RKC) , parameter :: TOL = sqrt(epsilon(0._RKC)) * 1000
299 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
300 : real(RKC) :: logx, output, output_ref, diff
301 : logical(LK) :: logMaxXPresent
302 :
303 7 : assertion = .true._LK
304 :
305 4008 : do i = 1_IK, 500_IK
306 4000 : logMinX = getUnifRand(-LOG_HUGE, LOG_HUGE)
307 4000 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
308 4000 : call runTestsWith(logMaxX = logMaxX)
309 4008 : call runTestsWith()
310 : end do
311 :
312 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 :
314 : contains
315 :
316 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
317 :
318 8000 : subroutine runTestsWith(logMaxX)
319 :
320 : use pm_val2str, only: getStr
321 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
322 : real(RKC), intent(in), optional :: logMaxX
323 :
324 8000 : logMaxXPresent = present(logMaxX)
325 8000 : alpha = -getUnifRand(0.5_RKC, +3._RKC)
326 8000 : if (logMaxXPresent) then
327 4000 : logCDFNF = getParetoLogCDFNF(alpha, logMinX, logMaxX)
328 : else
329 4000 : logCDFNF = 0._RKC ! getParetoLogCDFNF(alpha, logMinX)
330 : end if
331 :
332 :
333 8000 : logx = getUnifRand(logMinX, getOption(logMinX + 4._RKC, logMaxX))
334 8000 : if (logMaxXPresent) then
335 4000 : call setParetoLogCDF(output_ref, logx, alpha, logMinX, logCDFNF)
336 4000 : if (output_ref > 0._RKC) output_ref = 0._RKC ! ensure cdf cannot be larger than 1 due to roundoff error, otherwise algorithm checks will catch it.
337 : !write(*,*) output_ref, logx, alpha, logMinX
338 : #if getParetoLogQuan_ENABLED
339 2000 : output = getParetoLogQuan(output_ref, alpha, logMinX, logMaxX)
340 : #elif setParetoLogQuan_ENABLED
341 2000 : call setParetoLogQuan(output, output_ref, alpha, logMinX, logCDFNF)
342 : #else
343 : #error "Unrecognized interface."
344 : #endif
345 : else
346 4000 : call setParetoLogCDF(output_ref, logx, alpha, logMinX)
347 : !write(*,*) output_ref, logx, alpha, logMinX
348 : #if getParetoLogQuan_ENABLED
349 2000 : output = getParetoLogQuan(output_ref, alpha, logMinX)
350 : #elif setParetoLogQuan_ENABLED
351 2000 : call setParetoLogQuan(output, output_ref, alpha, logMinX)
352 : #endif
353 : end if
354 8000 : output_ref = logx
355 8000 : call report(__LINE__, logMaxX, logx)
356 :
357 8000 : end subroutine
358 :
359 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360 :
361 8000 : subroutine report(line, logMaxX, logx)
362 : integer, intent(in) :: line
363 : real(RKC), intent(in), optional :: logMaxX, logx
364 8000 : diff = abs(output - output_ref)
365 8000 : assertion = assertion .and. diff <= max(abs(output_ref) * TOL, TOL)
366 8000 : if (test%traceable .and. .not. assertion) then
367 : ! LCOV_EXCL_START
368 : write(test%disp%unit,"(*(g0,:,', '))")
369 : write(test%disp%unit,"(*(g0,:,', '))") "present(logx)", present(logx)
370 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMaxX)", present(logMaxX)
371 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, logx)", getOption(HUGE_RKC, logx)
372 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMaxX)", getOption(-HUGE_RKC, logMaxX)
373 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
374 : write(test%disp%unit,"(*(g0,:,', '))") "logMinX", logMinX
375 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
376 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
377 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
378 : write(test%disp%unit,"(*(g0,:,', '))") "output", output
379 : write(test%disp%unit,"(*(g0,:,', '))") "output_ref", output_ref
380 : write(test%disp%unit,"(*(g0,:,', '))")
381 : ! LCOV_EXCL_STOP
382 : end if
383 8000 : call test%assert(assertion, SK_"@getParetoLogQuan(): The Quantile of the Pareto distribution must be computed correctly.", int(line, IK))
384 8000 : end subroutine
385 :
386 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
387 :
388 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
389 : #elif getParetoLogRand_ENABLED || setParetoLogRand_ENABLED
390 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
391 :
392 : integer(IK) :: i
393 : integer(IK) , parameter :: NSIM = 1000_IK
394 : real(RKC) , parameter :: MEAN_REF = 0.5_RKC
395 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
396 : real(RKC) , parameter :: LOG_HUGE = log(HUGE_RKC)
397 : real(RKC) , parameter :: TOL = 0.1_RKC
398 : real(RKC) :: alpha, logCDFNF, logMinX, logMaxX
399 : real(RKC) :: LogRand(NSIM), CDF(NSIM), mean, diff
400 : logical(LK) :: logMaxXPresent
401 :
402 8 : assertion = .true._LK
403 :
404 24 : do i = 1_IK, 2_IK
405 16 : logMinX = getUnifRand(-LOG_HUGE, LOG_HUGE)
406 16 : logMaxX = getUnifRand(logMinX + 1._RKC, logMinX + 4._RKC)
407 16 : call runTestsWith(logMaxX = logMaxX)
408 24 : call runTestsWith()
409 : end do
410 :
411 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
412 :
413 : contains
414 :
415 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
416 :
417 32 : subroutine runTestsWith(logMaxX)
418 :
419 : use pm_val2str, only: getStr
420 : use pm_sampleMean, only: getMean
421 : use pm_distNegExp, only: getNegExpRand
422 : real(RKC), parameter :: SQRT_LOG_HUGE = sqrt(LOG_HUGE)
423 : real(RKC), intent(in), optional :: logMaxX
424 : integer(IK) :: j
425 :
426 32 : logMaxXPresent = present(logMaxX)
427 32 : alpha = -getUnifRand(0.5_RKC, +3._RKC)
428 32 : if (logMaxXPresent) then
429 16 : logCDFNF = getParetoLogCDFNF(alpha, logMinX, logMaxX)
430 : else
431 16 : logCDFNF = 0._RKC ! getParetoLogCDFNF(alpha, logMinX)
432 : end if
433 :
434 32032 : LogRand = getUnifRand(logMinX, getOption(logMinX + 4._RKC, logMaxX), s1 = size(LogRand, kind = IK))
435 32 : if (logMaxXPresent) then
436 16016 : do j = 1_IK, size(LogRand, kind = IK)
437 : #if getParetoLogRand_ENABLED
438 8000 : LogRand(j) = getParetoLogRand(alpha, logMinX, logMaxX) ! Truncated Pareto distribution.
439 : #elif setParetoLogRand_ENABLED
440 8000 : call setParetoLogRand(LogRand(j), getNegExpRand(1._RKC), alpha, logMinX, logCDFNF)
441 : #else
442 : #error "Unrecognized interface."
443 : #endif
444 16000 : call setParetoLogCDF(CDF(j), LogRand(j), alpha, logMinX, logCDFNF)
445 16016 : call report(__LINE__, logMaxX, LogRand(j))
446 : end do
447 : else
448 16016 : do j = 1_IK, size(LogRand, kind = IK)
449 : #if getParetoLogRand_ENABLED
450 8000 : LogRand(j) = getParetoLogRand(alpha, logMinX) ! Truncated Pareto distribution.
451 : #elif setParetoLogRand_ENABLED
452 8000 : call setParetoLogRand(LogRand(j), getNegExpRand(1._RKC), alpha, logMinX)
453 : #else
454 : #error "Unrecognized interface."
455 : #endif
456 16000 : call setParetoLogCDF(CDF(j), LogRand(j), alpha, logMinX)
457 16016 : call report(__LINE__, logMaxX, LogRand(j))
458 : end do
459 : end if
460 32032 : mean = getMean(exp(CDF))
461 32 : call report(__LINE__, logMaxX, mean = mean)
462 :
463 32 : end subroutine
464 :
465 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
466 :
467 32032 : subroutine report(line, logMaxX, logRand, mean)
468 : integer, intent(in) :: line
469 : real(RKC), intent(in), optional :: logMaxX, logRand, mean
470 32032 : if (present(mean) .and. present(logRand)) then
471 : error stop "Internal error occurred. Both `logRand` and `mean` input arguments cannot be present simultaneously." ! LCOV_EXCL_LINE
472 32032 : elseif (present(mean)) then
473 32 : diff = abs(mean - MEAN_REF)
474 32 : assertion = assertion .and. diff <= max(abs(MEAN_REF) * TOL, TOL)
475 32000 : elseif (present(logRand)) then
476 32000 : assertion = assertion .and. logical(logMinX <= logRand, LK) .and. logical(logRand <= getOption(logMaxX, huge(logMaxX)), LK)
477 : else
478 : error stop "Internal error occurred. Either `logRand` and `mean` input argument must be present." ! LCOV_EXCL_LINE
479 : end if
480 32032 : if (test%traceable .and. .not. assertion) then
481 : ! LCOV_EXCL_START
482 : write(test%disp%unit,"(*(g0,:,', '))")
483 : write(test%disp%unit,"(*(g0,:,', '))") "present(LogRand)", present(LogRand)
484 : write(test%disp%unit,"(*(g0,:,', '))") "present(logMaxX)", present(logMaxX)
485 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(HUGE_RKC, LogRand)", getOption(HUGE_RKC, LogRand)
486 : write(test%disp%unit,"(*(g0,:,', '))") "getOption(-HUGE_RKC, logMaxX)", getOption(-HUGE_RKC, logMaxX)
487 : write(test%disp%unit,"(*(g0,:,', '))") "logCDFNF", logCDFNF
488 : write(test%disp%unit,"(*(g0,:,', '))") "logMinX", logMinX
489 : write(test%disp%unit,"(*(g0,:,', '))") "alpha", alpha
490 : if (present(mean)) then
491 : write(test%disp%unit,"(*(g0,:,', '))") "diff", diff
492 : write(test%disp%unit,"(*(g0,:,', '))") "TOL", TOL
493 : write(test%disp%unit,"(*(g0,:,', '))") "mean", mean
494 : write(test%disp%unit,"(*(g0,:,', '))") "MEAN_REF", MEAN_REF
495 : write(test%disp%unit,"(*(g0,:,', '))")
496 : end if
497 : ! LCOV_EXCL_STOP
498 : end if
499 32032 : if (present(mean)) then
500 32 : call test%assert(assertion, SK_"The generated random value must be Pareto-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))
501 32000 : elseif (present(logRand)) then
502 32000 : call test%assert(assertion, SK_"The generated Pareto random value must be within the support of the distribution.", int(line, IK))
503 : end if
504 32032 : end subroutine
505 :
506 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
507 :
508 : #else
509 : !%%%%%%%%%%%%%%%%%%%%%%%%
510 : #error "Unrecognized interface."
511 : !%%%%%%%%%%%%%%%%%%%%%%%%
512 : #endif
|
