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 file contains procedure implementations of tests of [test_pm_sampleCCF](@ref test_pm_sampleCCF).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 5:03 PM, August 11, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 32 : character(:, SK), allocatable :: format
28 : ! Define the comparison precision and tolerance.
29 : real(TKC), parameter :: rtol = epsilon(1._TKC) * 5000
30 : #if CK_ENABLED
31 : #define TYPE_OF_SEQ complex(TKC)
32 : #define GET_CONJG(X)conjg(X)
33 : #define COERCE(X)X
34 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONES = (1._TKC, 1._TKC), TWOS = (2._TKC, 2._TKC)
35 : complex(TKC), parameter :: ctol = (rtol, rtol)
36 : #elif RK_ENABLED
37 : #define TYPE_OF_SEQ real(TKC)
38 : #define COERCE(X)real(X)
39 : #define GET_CONJG(X)X
40 : real(TKC), parameter :: ZERO = 0._TKC, ONES = 1._TKC, TWOS = 2._TKC
41 : real(TKC), parameter :: ctol = rtol
42 : #else
43 : #error "Unrecognized interface."
44 : #endif
45 : !%%%%%%%%%%%%%
46 : #if getCCF_ENABLED
47 : !%%%%%%%%%%%%%
48 :
49 : TYPE_OF_SEQ, allocatable :: f(:), g(:), ccf(:), ccf_ref(:), diff(:)
50 : integer(IK), allocatable :: lag(:)
51 : integer(IK) :: nf, ng, itry
52 : #if CK_ENABLED
53 4 : format = getFormat(deliml = '', subsep = SK_'', delimr = 'i', subcount = 2_IK, signed = .true._LK)
54 : #elif RK_ENABLED
55 4 : format = getFormat(mold = [0._TKC])
56 : #else
57 : #error "Unrecognized interface."
58 : #endif
59 8 : assertion = .true._LK
60 816 : do itry = 1, 100
61 :
62 : ! test FG pair interface.
63 :
64 : block
65 :
66 800 : nf = getUnifRand(2, 40)
67 800 : ng = getUnifRand(2, 40)
68 18303 : f = getUnifRand(ONES, TWOS * 3 / 2, nf)
69 18418 : g = getUnifRand(ONES, TWOS * 3 / 2, ng)
70 :
71 34321 : ccf = getCCF(f, g)
72 34321 : ccf_ref = getCCF_ref(f, g)
73 800 : call report(__LINE__)
74 :
75 67042 : ccf = getCCF(f, g, norm = zscore)
76 67042 : ccf_ref = getCCF_ref(f, g, norm = zscore)
77 800 : call report(__LINE__, norm = zscore)
78 :
79 67042 : ccf = getCCF(f, g, norm = nothing)
80 67042 : ccf_ref = getCCF_ref(f, g, norm = nothing)
81 800 : call report(__LINE__, norm = nothing)
82 :
83 67042 : ccf = getCCF(f, g, norm = stdscale)
84 67042 : ccf_ref = getCCF_ref(f, g, norm = stdscale)
85 800 : call report(__LINE__, norm = stdscale)
86 :
87 67042 : ccf = getCCF(f, g, norm = meanshift)
88 67042 : ccf_ref = getCCF_ref(f, g, norm = meanshift)
89 800 : call report(__LINE__, norm = meanshift)
90 :
91 : end block
92 :
93 : ! test FG pair interface with arbitrary `lag`.
94 :
95 8 : block
96 :
97 800 : nf = getUnifRand(2, 40)
98 800 : ng = getUnifRand(2, 40)
99 18869 : f = getUnifRand(ONES, TWOS * 3 / 2, nf)
100 18208 : g = getUnifRand(ONES, TWOS * 3 / 2, ng)
101 :
102 41600 : lag = getRange(1_IK, 50_IK)
103 :
104 41600 : ccf = getCCF(f, g, lag)
105 41600 : ccf_ref = getCCF_ref(f, g, lag)
106 800 : call report(__LINE__)
107 :
108 81600 : ccf = getCCF(f, g, lag, norm = zscore)
109 81600 : ccf_ref = getCCF_ref(f, g, lag, norm = zscore)
110 800 : call report(__LINE__, norm = zscore)
111 :
112 81600 : ccf = getCCF(f, g, lag, norm = nothing)
113 81600 : ccf_ref = getCCF_ref(f, g, lag, norm = nothing)
114 800 : call report(__LINE__, norm = nothing)
115 :
116 81600 : ccf = getCCF(f, g, lag, norm = stdscale)
117 81600 : ccf_ref = getCCF_ref(f, g, lag, norm = stdscale)
118 800 : call report(__LINE__, norm = stdscale)
119 :
120 81600 : ccf = getCCF(f, g, lag, norm = meanshift)
121 81600 : ccf_ref = getCCF_ref(f, g, lag, norm = meanshift)
122 800 : call report(__LINE__, norm = meanshift)
123 :
124 800 : lag = [integer(IK) ::]
125 1600 : ccf = getCCF(f, g, lag)
126 800 : ccf_ref = [TYPE_OF_SEQ ::]
127 800 : call report(__LINE__, lag)
128 :
129 : end block
130 :
131 : end do
132 :
133 : contains
134 :
135 8000 : function getCCF_ref(f, g, lag, norm) result(ccf)
136 : TYPE_OF_SEQ, intent(in), contiguous :: f(:), g(:)
137 : integer(IK), intent(in), contiguous, optional :: lag(:)
138 : TYPE_OF_SEQ, allocatable :: ff(:), gg(:), coef(:), ccf(:)
139 : class(*), intent(in), optional :: norm
140 : integer(IK), allocatable :: lag_def(:)
141 : integer(IK), allocatable :: factor(:)
142 16000 : class(*), allocatable :: norm_def
143 : real(TKC) :: dumf, dumg, normfac
144 : TYPE_OF_SEQ :: meanf, meang
145 : integer(IK) :: lagmin, lagmax
146 : integer(IK) :: lenccf, lenf, leng
147 : logical(LK) :: inf
148 8000 : lenf = size(f, 1, IK)
149 8000 : leng = size(g, 1, IK)
150 8000 : lagmax = leng - 1_IK
151 8000 : lagmin = 1_IK - lenf
152 8000 : if (present(lag)) then
153 4000 : if (size(lag, 1, IK) == 0_IK) then
154 0 : call setResized(ccf, 0_IK)
155 : return
156 : end if
157 4000 : lagmax = max(lagmax, lag(size(lag, 1, IK)))
158 4000 : lagmin = min(lagmin, lag(1))
159 208000 : lag_def = lag
160 : else
161 4776 : lag_def = getRange(lagmin, lagmax)
162 : end if
163 8000 : lenccf = lagmax - lagmin + 1_IK
164 8000 : call setResized(ff, lenccf)
165 8000 : call setResized(gg, lenccf)
166 8000 : call setResized(ccf, lenccf)
167 8000 : call setResized(coef, lenccf)
168 8000 : if (present(norm)) then
169 6400 : norm_def = norm
170 : else
171 1600 : norm_def = zscore
172 : end if
173 8000 : if (same_type_as(norm_def, meanshift) .or. same_type_as(norm_def, zscore)) then
174 4800 : meanf = getMean(f)
175 4800 : meang = getMean(g)
176 106716 : ff(1 : lenf) = f - meanf
177 105078 : gg(1 : leng) = g - meang
178 : else
179 71144 : ff(1 : lenf) = f
180 70052 : gg(1 : leng) = g
181 : end if
182 288090 : ff(lenf + 1 :) = ZERO
183 290820 : gg(leng + 1 :) = ZERO
184 8000 : normfac = 1._TKC / lenccf
185 8000 : if (same_type_as(norm_def, stdscale) .or. same_type_as(norm_def, zscore)) then
186 106716 : dumf = real(dot_product(ff(1 : lenf), ff(1 : lenf)), TKC)
187 105078 : dumg = real(dot_product(gg(1 : leng), gg(1 : leng)), TKC)
188 4800 : normfac = normfac / sqrt(dumf * dumg)
189 : end if
190 33550 : factor = getFactorFFT(ff, coef)
191 8000 : call setCCF(factor, coef, ff, gg, ccf, inf)
192 8000 : if (inf) then
193 215660 : ccf = ff * normfac
194 : else
195 250290 : ccf = gg * normfac
196 : end if
197 915900 : ccf = cshift(ccf, lagmin)
198 1106815 : ccf = ccf(lag_def - lagmin + 1)
199 8000 : end function
200 :
201 8800 : subroutine report(line, lag, norm)
202 : integer, intent(in) :: line
203 : class(*), intent(in), optional :: norm
204 : integer(IK), intent(in), optional :: lag(:)
205 8800 : if (0_IK < size(ccf, 1, IK)) then
206 560110 : diff = abs(ccf - ccf_ref)
207 371605 : assertion = assertion .and. all(diff < ctol)
208 : else
209 800 : assertion = assertion .and. size(ccf_ref, 1, IK) == 0_IK
210 : end if
211 8800 : if (test%traceable .and. .not. assertion) then
212 : ! LCOV_EXCL_START
213 : call test%disp%skip()
214 : call test%disp%show("[nf, ng, size(ccf, 1, IK)]")
215 : call test%disp%show( [nf, ng, size(ccf, 1, IK)] )
216 : call test%disp%show("f")
217 : call test%disp%show( f , format = format )
218 : call test%disp%show("g")
219 : call test%disp%show( g , format = format )
220 : call test%disp%show("ccf")
221 : call test%disp%show( ccf , format = format )
222 : call test%disp%show("ccf_ref")
223 : call test%disp%show( ccf_ref , format = format )
224 : call test%disp%show("diff")
225 : call test%disp%show( diff , format = format )
226 : call test%disp%show("present(lag)")
227 : call test%disp%show( present(lag) )
228 : if (present(lag)) then
229 : call test%disp%show("lag")
230 : call test%disp%show( lag )
231 : end if
232 : call test%disp%show("present(norm)")
233 : call test%disp%show( present(norm) )
234 : if (present(norm)) then
235 : call test%disp%show("[same_type_as(norm, meanshift), same_type_as(norm, stdscale), same_type_as(norm, zscore), same_type_as(norm, nothing)]")
236 : call test%disp%show( [same_type_as(norm, meanshift), same_type_as(norm, stdscale), same_type_as(norm, zscore), same_type_as(norm, nothing)] )
237 : end if
238 : call test%disp%skip()
239 : ! LCOV_EXCL_STOP
240 : end if
241 8800 : call test%assert(assertion, SK_"The `ccf` must be correctly computed for the specified sequences.", int(line, IK))
242 8800 : end subroutine
243 :
244 : !%%%%%%%%%%%%%
245 : #elif setCCF_ENABLED
246 : !%%%%%%%%%%%%%
247 :
248 : TYPE_OF_SEQ, allocatable :: f(:), g(:), ff(:), gg(:), coef(:), ccf(:), ccf_ref(:), diff(:)
249 : integer(IK), allocatable :: factor(:)
250 : integer(IK) :: nf, ng, nc, itry
251 : logical(LK) :: inf
252 : #if CK_ENABLED
253 4 : format = getFormat(deliml = '', subsep = SK_'', delimr = 'i', subcount = 2_IK, signed = .true._LK)
254 : #elif RK_ENABLED
255 4 : format = getFormat(mold = [0._TKC])
256 : #else
257 : #error "Unrecognized interface."
258 : #endif
259 8 : assertion = .true._LK
260 176 : do itry = 2, 42, 2
261 :
262 : ! test FG pair interface.
263 :
264 8 : block
265 168 : nf = itry
266 168 : ng = itry + 1
267 168 : nc = nf + ng - 1_IK
268 :
269 168 : call setResized(f, nc)
270 168 : call setResized(g, nc)
271 3864 : f(1:nf) = getUnifRand(ONES, TWOS * 3 / 2, nf)
272 4032 : g(1:ng) = getUnifRand(ONES, TWOS * 3 / 2, ng)
273 3864 : f(nf + 1:) = ZERO
274 3696 : g(ng + 1:) = ZERO
275 7728 : ccf_ref = getCCF_ref(f, g)
276 :
277 7896 : ff = f
278 7896 : gg = g
279 168 : call setResized(ccf, nc)
280 168 : call setResized(coef, nc)
281 760 : factor = getFactorFFT(f, coef)
282 168 : call setCCF(factor, coef, ff, gg, ccf, inf)
283 7728 : ccf = merge(ff / nc, gg / nc, inf)
284 168 : call report(__LINE__)
285 : end block
286 :
287 : end do
288 :
289 : contains
290 :
291 168 : function getCCF_ref(f, g) result(ccf)
292 : TYPE_OF_SEQ, intent(in) :: f(:), g(:)
293 : complex(TKC), allocatable :: ff(:), gg(:)
294 : TYPE_OF_SEQ, allocatable :: ccf(:)
295 26544 : ff = f; gg = g; ccf = COERCE(getFFTI(getFFTF(ff) * conjg(getFFTF(gg))))
296 168 : end function
297 :
298 168 : subroutine report(line)
299 : integer, intent(in) :: line
300 11424 : diff = abs(ccf - ccf_ref)
301 7560 : assertion = assertion .and. all(diff < ctol)
302 168 : if (test%traceable .and. .not. assertion) then
303 : ! LCOV_EXCL_START
304 : call test%disp%skip()
305 : call test%disp%show("[nf, ng, nc]")
306 : call test%disp%show( [nf, ng, nc] )
307 : call test%disp%show("f")
308 : call test%disp%show( f , format = format )
309 : call test%disp%show("g")
310 : call test%disp%show( g , format = format )
311 : call test%disp%show("ccf")
312 : call test%disp%show( ccf , format = format )
313 : call test%disp%show("ccf_ref")
314 : call test%disp%show( ccf_ref , format = format )
315 : call test%disp%show("diff")
316 : call test%disp%show( diff , format = format )
317 : call test%disp%skip()
318 : ! LCOV_EXCL_STOP
319 : end if
320 168 : call test%assert(assertion, SK_"The `ccf` must be correctly computed for the specified sequences.", int(line, IK))
321 168 : end subroutine
322 :
323 : !%%%%%%%%%%%%%
324 : #elif getACF_ENABLED
325 : !%%%%%%%%%%%%%
326 :
327 : TYPE_OF_SEQ, allocatable :: f(:), g(:), acf(:), acf_ref(:), diff(:)
328 : integer(IK), allocatable :: lag(:)
329 : integer(IK) :: nf, itry
330 : #if CK_ENABLED
331 4 : format = getFormat(deliml = '', subsep = SK_'', delimr = 'i', subcount = 2_IK, signed = .true._LK)
332 : #elif RK_ENABLED
333 4 : format = getFormat(mold = [0._TKC])
334 : #else
335 : #error "Unrecognized interface."
336 : #endif
337 8 : assertion = .true._LK
338 :
339 816 : do itry = 1, 100
340 :
341 : ! test FG pair interface.
342 :
343 800 : nf = getUnifRand(2, 40)
344 18433 : f = getUnifRand(ONES, TWOS * 3 / 2, nf)
345 :
346 800 : call runTestWith(__LINE__)
347 800 : call runTestWith(__LINE__, norm = zscore)
348 800 : call runTestWith(__LINE__, norm = stdscale)
349 800 : call runTestWith(__LINE__, norm = meanshift)
350 :
351 1600 : lag = getRange(1_IK, 50_IK)
352 800 : call runTestWith(__LINE__, lag)
353 800 : call runTestWith(__LINE__, lag, norm = zscore)
354 800 : call runTestWith(__LINE__, lag, norm = stdscale)
355 800 : call runTestWith(__LINE__, lag, norm = meanshift)
356 :
357 800 : lag = [integer(IK) ::]
358 800 : call runTestWith(__LINE__, lag)
359 800 : call runTestWith(__LINE__, lag, norm = zscore)
360 800 : call runTestWith(__LINE__, lag, norm = stdscale)
361 808 : call runTestWith(__LINE__, lag, norm = meanshift)
362 :
363 : end do
364 :
365 : contains
366 :
367 9600 : subroutine runTestWith(line, lag, norm)
368 : integer, intent(in) :: line
369 : class(*), intent(in), optional :: norm
370 : integer(IK), intent(in), contiguous, optional :: lag(:)
371 : integer(IK), allocatable :: lag_def(:)
372 249732 : acf = getACF(f, lag, norm)
373 221974 : g = f
374 9600 : if (present(lag)) then
375 172800 : acf_ref = getCCF(f, g, lag, norm)
376 : else
377 3978 : lag_def = getRange(0_IK, size(f, 1, IK) - 1_IK)
378 73732 : acf_ref = getCCF(f, g, lag_def, norm)
379 : end if
380 12800 : call report(line, lag, norm)
381 9600 : end subroutine
382 :
383 9600 : subroutine report(line, lag, norm)
384 : integer, intent(in) :: line
385 : class(*), intent(in), optional :: norm
386 : integer(IK), intent(in), optional :: lag(:)
387 9600 : if (0_IK < size(acf, 1, IK)) then
388 353512 : diff = abs(acf - acf_ref)
389 233732 : assertion = assertion .and. all(diff < ctol)
390 : else
391 3200 : assertion = assertion .and. size(acf_ref, 1, IK) == 0_IK
392 : end if
393 9600 : if (test%traceable .and. .not. assertion) then
394 : ! LCOV_EXCL_START
395 : call test%disp%skip()
396 : call test%disp%show("[nf, size(acf, 1, IK)]")
397 : call test%disp%show( [nf, size(acf, 1, IK)] )
398 : call test%disp%show("f")
399 : call test%disp%show( f , format = format )
400 : call test%disp%show("acf")
401 : call test%disp%show( acf , format = format )
402 : call test%disp%show("acf_ref")
403 : call test%disp%show( acf_ref , format = format )
404 : call test%disp%show("diff")
405 : call test%disp%show( diff , format = format )
406 : call test%disp%show("present(lag)")
407 : call test%disp%show( present(lag) )
408 : if (present(lag)) then
409 : call test%disp%show("lag")
410 : call test%disp%show( lag )
411 : end if
412 : call test%disp%show("present(norm)")
413 : call test%disp%show( present(norm) )
414 : if (present(norm)) then
415 : call test%disp%show("[same_type_as(norm, meanshift), same_type_as(norm, stdscale), same_type_as(norm, zscore), same_type_as(norm, nothing)]")
416 : call test%disp%show( [same_type_as(norm, meanshift), same_type_as(norm, stdscale), same_type_as(norm, zscore), same_type_as(norm, nothing)] )
417 : end if
418 : call test%disp%skip()
419 : ! LCOV_EXCL_STOP
420 : end if
421 9600 : call test%assert(assertion, SK_"The `acf` must be correctly computed for the specified sequences.", int(line, IK))
422 9600 : end subroutine
423 :
424 : !%%%%%%%%%%%%%
425 : #elif setACF_ENABLED
426 : !%%%%%%%%%%%%%
427 :
428 : TYPE_OF_SEQ, allocatable :: seq(:), f(:), g(:), coef(:), acf(:), acf_ref(:), diff(:)
429 : integer(IK), allocatable :: factor(:)
430 : integer(IK) :: nf, nc, itry
431 : logical(LK) :: inf
432 : #if CK_ENABLED
433 4 : format = getFormat(deliml = '', subsep = SK_'', delimr = 'i', subcount = 2_IK, signed = .true._LK)
434 : #elif RK_ENABLED
435 4 : format = getFormat(mold = [0._TKC])
436 : #else
437 : #error "Unrecognized interface."
438 : #endif
439 8 : assertion = .true._LK
440 176 : do itry = 2, 42, 2
441 :
442 : ! test FG pair interface.
443 :
444 8 : block
445 168 : nf = itry
446 168 : nc = nf * 2_IK - 1_IK
447 168 : call setResized(seq, nc)
448 3864 : seq(1:nf) = getUnifRand(ONES, TWOS * 3 / 2, nf)
449 3696 : seq(nf + 1:) = ZERO
450 7728 : f = seq
451 7728 : g = seq
452 :
453 168 : call setResized(acf, nc)
454 168 : call setResized(coef, nc)
455 168 : call setResized(acf_ref, nc)
456 592 : factor = getFactorFFT(f, coef)
457 168 : call setCCF(factor, coef, f, g, acf_ref, inf)
458 7560 : acf_ref = merge(f, g, inf)
459 7560 : f = seq
460 168 : call setACF(factor, coef, f, acf, inf)
461 7560 : acf = merge(f, acf, inf)
462 168 : call report(__LINE__)
463 : end block
464 :
465 : end do
466 :
467 : contains
468 :
469 168 : subroutine report(line)
470 : integer, intent(in) :: line
471 11172 : diff = abs(acf - acf_ref)
472 7392 : assertion = assertion .and. all(diff < ctol)
473 168 : if (test%traceable .and. .not. assertion) then
474 : ! LCOV_EXCL_START
475 : call test%disp%skip()
476 : call test%disp%show("[nf, nc]")
477 : call test%disp%show( [nf, nc] )
478 : call test%disp%show("seq")
479 : call test%disp%show( seq , format = format )
480 : call test%disp%show("acf")
481 : call test%disp%show( acf , format = format )
482 : call test%disp%show("acf_ref")
483 : call test%disp%show( acf_ref , format = format )
484 : call test%disp%show("diff")
485 : call test%disp%show( diff , format = format )
486 : call test%disp%skip()
487 : ! LCOV_EXCL_STOP
488 : end if
489 168 : call test%assert(assertion, SK_"The `acf` must be correctly computed for the specified sequence.", int(line, IK))
490 168 : end subroutine
491 :
492 : #else
493 : !%%%%%%%%%%%%%%%%%%%%%%%%
494 : #error "Unrecognized interface."
495 : !%%%%%%%%%%%%%%%%%%%%%%%%
496 : #endif
497 : #undef TYPE_OF_SEQ
498 : #undef GET_CONJG
499 : #undef COERCE
|
