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 [pm_sampleCor](@ref pm_sampleCor).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 5:03 PM, August 11, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Define the comparison precision and tolerance.
28 : #if CK_ENABLED && (getCor_ENABLED || setCor_ENABLED)
29 : real(TKC), parameter :: rtol = epsilon(1._TKC) * 100
30 : complex(TKC), parameter :: ctol = (rtol, rtol)
31 : #define GET_CONJG(X)conjg(X)
32 : #elif RK_ENABLED && (getCor_ENABLED || setCor_ENABLED)
33 : real(TKC), parameter :: rtol = epsilon(1._TKC) * 100
34 : real(TKC), parameter :: ctol = rtol
35 : #define GET_CONJG(X)X
36 : #elif getCor_ENABLED || setCor_ENABLED
37 : #error "Unrecognized interface."
38 : #else
39 : real(RK), parameter :: rtol = epsilon(1._RK) * 100
40 : #endif
41 : ! Define the sample type.
42 : #if SK_ENABLED && D0_ENABLED
43 : #define GET_SHAPE len
44 : #define TYPE_OF_SAMPLE character(:,TKC)
45 : character(1,TKC), parameter :: lb = TKC_"a", ub = TKC_"z"
46 : #else
47 : #define GET_SHAPE shape
48 : #if SK_ENABLED
49 : #define TYPE_OF_SAMPLE character(2,TKC)
50 : TYPE_OF_SAMPLE, parameter :: lb = TKC_"aa", ub = TKC_"zz"
51 : #elif IK_ENABLED
52 : #define TYPE_OF_SAMPLE integer(TKC)
53 : TYPE_OF_SAMPLE, parameter :: lb = -9_TKC, ub = +9_TKC
54 : #elif CK_ENABLED
55 : #define TYPE_OF_SAMPLE complex(TKC)
56 : TYPE_OF_SAMPLE, parameter :: lb = (2._TKC, -3._TKC), ub = (3._TKC, -2._TKC), ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC), ONES = (1._TKC, 1._TKC), TWOS = ONES + ONES
57 : #elif RK_ENABLED
58 : #define TYPE_OF_SAMPLE real(TKC)
59 : TYPE_OF_SAMPLE, parameter :: lb = 2._TKC, ub = 3._TKC, ZERO = 0._TKC, ONE = 1._TKC, ONES = 1._TKC, TWOS = ONES + ONES
60 : #elif PSSK_ENABLED
61 : #define TYPE_OF_SAMPLE type(css_pdt(TKC))
62 : character(1,TKC), parameter :: lb = TKC_"a", ub = TKC_"z"
63 : #elif BSSK_ENABLED
64 : #define TYPE_OF_SAMPLE type(css_type)
65 : character(1,TKC), parameter :: lb = TKC_"a", ub = TKC_"z"
66 : #elif !setCordance_ENABLED
67 : #error "Unrecognized interface."
68 : #endif
69 : #endif
70 : !%%%%%%%%%%%%%
71 : #if getCor_ENABLED
72 : !%%%%%%%%%%%%%
73 :
74 : real(TKC) :: rweisum
75 : integer(IK) :: iweisum
76 : real(TKC), allocatable :: rweight(:)
77 : integer(IK), allocatable :: iweight(:)
78 : integer(IK) :: itry, nsam, ndim, dim
79 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:)
80 : TYPE_OF_SAMPLE, allocatable :: cor(:,:), cov(:,:), covupp(:,:), covlow(:,:), cor_ref(:,:), cdiff(:,:)
81 : real(TKC), allocatable :: std(:), stdinv(:)
82 7 : assertion = .true._LK
83 :
84 408 : do itry = 1, 50
85 :
86 : ! Test CFC: The strategy is to construct a cot mat, convert to cov mat using `setCov`, and recover the original cor mat using `setCor`.
87 :
88 : block
89 :
90 400 : ndim = getUnifRand(1_IK, 5_IK)
91 6210 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
92 400 : call setMatCopy(cor_ref, rdpack, cor_ref, rdpack, upp, transHerm)
93 400 : call setMatInit(cor_ref, dia, ONE)
94 1969 : std = getUnifRand(1._TKC, 2._TKC, ndim)
95 1969 : stdinv = 1._TKC / std
96 11620 : cov = getCov(cor_ref, uppDia, std)
97 6862 : covupp = cov; call setMatInit(covupp, low, ZERO)
98 6862 : covlow = cov; call setMatInit(covlow, upp, ZERO)
99 :
100 11620 : cor = getCor(covupp, uppDia)
101 400 : call reportCFC(__LINE__)
102 :
103 11620 : cor = getCor(covlow, lowDia)
104 400 : call reportCFC(__LINE__)
105 :
106 400 : call setMatInit(covupp, dia, ZERO)
107 400 : call setMatInit(covlow, dia, ZERO)
108 :
109 11620 : cor = getCor(covupp, upp, stdinv)
110 400 : call reportCFC(__LINE__, stdinv)
111 :
112 11620 : cor = getCor(covlow, low, stdinv)
113 400 : call reportCFC(__LINE__, stdinv)
114 :
115 : end block
116 :
117 : ! Test Prs: The strategy is to construct a sample, compute its correlation and compare it against the reference computed here.
118 :
119 400 : nsam = getUnifRand(2_IK, 7_IK)
120 :
121 : ! test sample interface.
122 :
123 : block
124 :
125 : type(uppDia_type), parameter :: subsetr = uppDia_type()
126 400 : dim = getChoice([1, 2])
127 400 : ndim = getUnifRand(1_IK, 5_IK)
128 400 : if (dim == 2) then
129 4217 : sample = getUnifRand(ONES, TWOS, ndim, nsam)
130 : else
131 3611 : sample = getUnifRand(ONES, TWOS, nsam, ndim)
132 : end if
133 2625 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
134 2625 : iweight = getUnifRand(1_IK, 9_IK, nsam)
135 2225 : iweisum = sum(iweight)
136 2225 : rweisum = sum(rweight)
137 :
138 6244 : cor_ref = getFilled(ZERO, ndim, ndim)
139 6244 : cor = getFilled(ZERO, ndim, ndim)
140 :
141 : ! integer weighted
142 :
143 1972 : mean = getMean(sample, dim, iweight)
144 :
145 6244 : cor = getCor(sample, dim, iweight)
146 400 : call setCor(cor_ref, subsetr, mean, sample, dim, iweight, iweisum)
147 400 : call setMatCopy(cor_ref, rdpack, cor_ref, rdpack, subsetr, transHerm)
148 400 : call report(__LINE__, SK_"integer-weighted")
149 :
150 : ! real weighted
151 :
152 1972 : mean = getMean(sample, dim, rweight)
153 :
154 6244 : cor = getCor(sample, dim, rweight)
155 400 : call setCor(cor_ref, subsetr, mean, sample, dim, rweight, rweisum)
156 400 : call setMatCopy(cor_ref, rdpack, cor_ref, rdpack, subsetr, transHerm)
157 400 : call report(__LINE__, SK_"real-weighted")
158 :
159 : ! unweighted
160 :
161 1972 : mean = getMean(sample, dim)
162 :
163 6244 : cor = getCor(sample, dim)
164 400 : call setCor(cor_ref, subsetr, mean, sample, dim)
165 400 : call setMatCopy(cor_ref, rdpack, cor_ref, rdpack, subsetr, transHerm)
166 400 : call report(__LINE__, SK_"unweighted")
167 :
168 : end block
169 :
170 : ! test sample XY interface.
171 :
172 8 : block
173 :
174 400 : dim = 1
175 400 : ndim = 2
176 5250 : sample = getUnifRand(ONES, TWOS, nsam, ndim)
177 2625 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
178 2625 : iweight = getUnifRand(1_IK, 9_IK, nsam)
179 2225 : iweisum = sum(iweight)
180 2225 : rweisum = sum(rweight)
181 :
182 3200 : cor_ref = getFilled(ZERO, ndim, ndim)
183 3200 : cor = getFilled(ZERO, ndim, ndim)
184 :
185 : ! integer weighted
186 :
187 1600 : mean = getMean(sample, dim, iweight)
188 :
189 400 : cor(1,2) = getCor(sample(:,1), sample(:,2), iweight)
190 400 : call setCor(cor_ref(1,2), mean, sample(:,1), sample(:,2), iweight, iweisum)
191 400 : call report(__LINE__, SK_"integer-weighted")
192 :
193 : ! real weighted
194 :
195 1600 : mean = getMean(sample, dim, rweight)
196 :
197 400 : cor(1,2) = getCor(sample(:,1), sample(:,2), rweight)
198 400 : call setCor(cor_ref(1,2), mean, sample(:,1), sample(:,2), rweight, rweisum)
199 400 : call report(__LINE__, SK_"real-weighted")
200 :
201 : ! unweighted
202 :
203 1600 : mean = getMean(sample, dim)
204 :
205 400 : cor(1,2) = getCor(sample(:,1), sample(:,2))
206 400 : call setCor(cor_ref(1,2), mean, sample(:,1), sample(:,2))
207 400 : call report(__LINE__, SK_"unweighted")
208 :
209 : end block
210 :
211 : end do
212 :
213 : contains
214 :
215 1600 : subroutine reportCFC(line, stdinv)
216 : real(TKC), intent(in), optional :: stdinv(:)
217 : integer, intent(in) :: line
218 36176 : cdiff = abs(cor - cor_ref)
219 23240 : assertion = assertion .and. all(cdiff < ctol)
220 1600 : if (test%traceable .and. .not. assertion) then
221 : ! LCOV_EXCL_START
222 : call test%disp%skip()
223 : call test%disp%show("ndim")
224 : call test%disp%show( ndim )
225 : call test%disp%show("present(stdinv)")
226 : call test%disp%show( present(stdinv) )
227 : if (present(stdinv)) then
228 : call test%disp%show("stdinv")
229 : call test%disp%show( stdinv )
230 : end if
231 : call test%disp%show("cov")
232 : call test%disp%show( cov )
233 : call test%disp%show("cor")
234 : call test%disp%show( cor )
235 : call test%disp%show("cor_ref")
236 : call test%disp%show( cor_ref )
237 : call test%disp%show("cdiff")
238 : call test%disp%show( cdiff )
239 : call test%disp%skip()
240 : ! LCOV_EXCL_STOP
241 : end if
242 1600 : call test%assert(assertion, SK_"The `cor` must be correctly computed from the specified `cor` and `stdinv`.", int(line, IK))
243 1600 : end subroutine
244 :
245 2400 : subroutine report(line, this)
246 : integer, intent(in) :: line
247 : character(*, SK), intent(in) :: this
248 39648 : cdiff = abs(cor - cor_ref)
249 25932 : assertion = assertion .and. all(cdiff < ctol)
250 2400 : if (test%traceable .and. .not. assertion) then
251 : ! LCOV_EXCL_START
252 : call test%disp%skip()
253 : call test%disp%show("[ndim, nsam, dim, shape(sample, IK)]")
254 : call test%disp%show( [ndim, nsam, dim, shape(sample, IK)] )
255 : call test%disp%show("sample")
256 : call test%disp%show( sample )
257 : call test%disp%show("iweight")
258 : call test%disp%show( iweight )
259 : call test%disp%show("rweight")
260 : call test%disp%show( rweight )
261 : call test%disp%show("cor_ref")
262 : call test%disp%show( cor_ref )
263 : call test%disp%show("cor")
264 : call test%disp%show( cor )
265 : call test%disp%show("cdiff")
266 : call test%disp%show( cdiff )
267 : call test%disp%skip()
268 : ! LCOV_EXCL_STOP
269 : end if
270 2400 : call test%assert(assertion, SK_"The `cor` must be computed correctly for "//this//SK_" sample.", int(line, IK))
271 2400 : end subroutine
272 :
273 : !%%%%%%%%%%%%%
274 : #elif setCor_ENABLED
275 : !%%%%%%%%%%%%%
276 :
277 : real(TKC) :: rweisum
278 : integer(IK) :: iweisum
279 : real(TKC), allocatable :: rweight(:)
280 : integer(IK), allocatable :: iweight(:)
281 : integer(IK) :: itry, nsam, ndim, dim
282 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:)
283 : TYPE_OF_SAMPLE, allocatable :: cor(:,:), cov(:,:), cor_ref(:,:), cdiff(:,:)
284 : real(TKC), allocatable :: std(:), stdinv(:)
285 8 : assertion = .true._LK
286 :
287 408 : do itry = 1, 50
288 :
289 : ! Test CFC: The strategy is to construct a cot mat, convert to cov mat using `setCov`, and recover the original cor mat using `setCor`.
290 :
291 : block
292 :
293 : block
294 :
295 : type(low_type) :: subsetr
296 : type(low_type) :: subsetv
297 400 : ndim = getUnifRand(1_IK, 5_IK)
298 6600 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
299 400 : call setMatInit(cor_ref, getSubComp(subsetr), ZERO, ZERO)
300 2032 : std = getUnifRand(1._TKC, 2._TKC, ndim)
301 2032 : stdinv = 1._TKC / std
302 6600 : cov = getFilled(ZERO, ndim, ndim)
303 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
304 :
305 6600 : cor = getFilled(ZERO, ndim, ndim)
306 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
307 400 : call reportCFC(__LINE__)
308 :
309 400 : call setMatInit(cov, dia, ZERO)
310 6600 : cor = getFilled(ZERO, ndim, ndim)
311 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
312 400 : call reportCFC(__LINE__, stdinv)
313 :
314 : end block
315 :
316 : block
317 :
318 : type(low_type) :: subsetr
319 : type(upp_type) :: subsetv
320 400 : ndim = getUnifRand(1_IK, 5_IK)
321 6458 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
322 400 : call setMatInit(cor_ref, getSubComp(subsetr), ZERO, ZERO)
323 2011 : std = getUnifRand(1._TKC, 2._TKC, ndim)
324 2011 : stdinv = 1._TKC / std
325 6458 : cov = getFilled(ZERO, ndim, ndim)
326 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
327 :
328 6458 : cor = getFilled(ZERO, ndim, ndim)
329 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
330 400 : call reportCFC(__LINE__)
331 :
332 400 : call setMatInit(cov, dia, ZERO)
333 6458 : cor = getFilled(ZERO, ndim, ndim)
334 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
335 400 : call reportCFC(__LINE__, stdinv)
336 :
337 : end block
338 :
339 : block
340 :
341 : type(upp_type) :: subsetr
342 : type(low_type) :: subsetv
343 400 : ndim = getUnifRand(1_IK, 5_IK)
344 6350 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
345 400 : call setMatInit(cor_ref, getSubComp(subsetr), ZERO, ZERO)
346 1992 : std = getUnifRand(1._TKC, 2._TKC, ndim)
347 1992 : stdinv = 1._TKC / std
348 6350 : cov = getFilled(ZERO, ndim, ndim)
349 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
350 :
351 6350 : cor = getFilled(ZERO, ndim, ndim)
352 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
353 400 : call reportCFC(__LINE__)
354 :
355 400 : call setMatInit(cov, dia, ZERO)
356 6350 : cor = getFilled(ZERO, ndim, ndim)
357 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
358 400 : call reportCFC(__LINE__, stdinv)
359 :
360 : end block
361 :
362 : block
363 :
364 : type(upp_type) :: subsetr
365 : type(upp_type) :: subsetv
366 400 : ndim = getUnifRand(1_IK, 5_IK)
367 6476 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
368 400 : call setMatInit(cor_ref, getSubComp(subsetr), ZERO, ZERO)
369 2013 : std = getUnifRand(1._TKC, 2._TKC, ndim)
370 2013 : stdinv = 1._TKC / std
371 6476 : cov = getFilled(ZERO, ndim, ndim)
372 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
373 :
374 6476 : cor = getFilled(ZERO, ndim, ndim)
375 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
376 400 : call reportCFC(__LINE__)
377 :
378 400 : call setMatInit(cov, dia, ZERO)
379 6476 : cor = getFilled(ZERO, ndim, ndim)
380 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
381 400 : call reportCFC(__LINE__, stdinv)
382 :
383 : end block
384 :
385 : block
386 :
387 : type(lowDia_type) :: subsetr
388 : type(low_type) :: subsetv
389 400 : ndim = getUnifRand(1_IK, 5_IK)
390 6104 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
391 400 : call setMatInit(cor_ref, getSubSymm(subsetr), ZERO, ONE)
392 1955 : std = getUnifRand(1._TKC, 2._TKC, ndim)
393 1955 : stdinv = 1._TKC / std
394 6104 : cov = getFilled(ZERO, ndim, ndim)
395 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
396 :
397 6104 : cor = getFilled(ZERO, ndim, ndim)
398 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
399 400 : call reportCFC(__LINE__)
400 :
401 400 : call setMatInit(cov, dia, ZERO)
402 6104 : cor = getFilled(ZERO, ndim, ndim)
403 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
404 400 : call reportCFC(__LINE__, stdinv)
405 :
406 : end block
407 :
408 : block
409 :
410 : type(lowDia_type) :: subsetr
411 : type(upp_type) :: subsetv
412 400 : ndim = getUnifRand(1_IK, 5_IK)
413 6028 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
414 400 : call setMatInit(cor_ref, getSubSymm(subsetr), ZERO, ONE)
415 1948 : std = getUnifRand(1._TKC, 2._TKC, ndim)
416 1948 : stdinv = 1._TKC / std
417 6028 : cov = getFilled(ZERO, ndim, ndim)
418 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
419 :
420 6028 : cor = getFilled(ZERO, ndim, ndim)
421 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
422 400 : call reportCFC(__LINE__)
423 :
424 400 : call setMatInit(cov, dia, ZERO)
425 6028 : cor = getFilled(ZERO, ndim, ndim)
426 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
427 400 : call reportCFC(__LINE__, stdinv)
428 :
429 : end block
430 :
431 : block
432 :
433 : type(uppDia_type) :: subsetr
434 : type(low_type) :: subsetv
435 400 : ndim = getUnifRand(1_IK, 5_IK)
436 6546 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
437 400 : call setMatInit(cor_ref, getSubSymm(subsetr), ZERO, ONE)
438 2014 : std = getUnifRand(1._TKC, 2._TKC, ndim)
439 2014 : stdinv = 1._TKC / std
440 6546 : cov = getFilled(ZERO, ndim, ndim)
441 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
442 :
443 6546 : cor = getFilled(ZERO, ndim, ndim)
444 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
445 400 : call reportCFC(__LINE__)
446 :
447 400 : call setMatInit(cov, dia, ZERO)
448 6546 : cor = getFilled(ZERO, ndim, ndim)
449 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
450 400 : call reportCFC(__LINE__, stdinv)
451 :
452 : end block
453 :
454 : block
455 :
456 : type(uppDia_type) :: subsetr
457 : type(upp_type) :: subsetv
458 400 : ndim = getUnifRand(1_IK, 5_IK)
459 6422 : cor_ref = getUnifRand(-ONES, ONES, ndim, ndim)
460 400 : call setMatInit(cor_ref, getSubSymm(subsetr), ZERO, ONE)
461 2002 : std = getUnifRand(1._TKC, 2._TKC, ndim)
462 2002 : stdinv = 1._TKC / std
463 6422 : cov = getFilled(ZERO, ndim, ndim)
464 400 : call setCov(cov, getSubUnion(subsetv, dia), cor_ref, getSubUnion(subsetr, dia), std)
465 :
466 6422 : cor = getFilled(ZERO, ndim, ndim)
467 400 : call setCor(cor, subsetr, cov, getSubUnion(subsetv, dia))
468 400 : call reportCFC(__LINE__)
469 :
470 400 : call setMatInit(cov, dia, ZERO)
471 6422 : cor = getFilled(ZERO, ndim, ndim)
472 400 : call setCor(cor, subsetr, cov, subsetv, stdinv)
473 400 : call reportCFC(__LINE__, stdinv)
474 :
475 : end block
476 : end block
477 :
478 : ! Test Prs: The strategy is to construct a sample, compute its correlation and compare it against the reference computed here.
479 :
480 400 : nsam = getUnifRand(2_IK, 7_IK)
481 :
482 : ! test sample lowDia interface.
483 :
484 : block
485 :
486 : type(lowDia_type), parameter :: subsetr = lowDia_type()
487 400 : dim = getChoice([1, 2])
488 400 : ndim = getUnifRand(1_IK, 5_IK)
489 400 : if (dim == 2) then
490 3975 : sample = getUnifRand(ONES, TWOS, ndim, nsam)
491 : else
492 3641 : sample = getUnifRand(ONES, TWOS, nsam, ndim)
493 : end if
494 2582 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
495 2582 : iweight = getUnifRand(1_IK, 9_IK, nsam)
496 2182 : iweisum = sum(iweight)
497 2182 : rweisum = sum(rweight)
498 :
499 : ! integer weighted
500 :
501 1995 : mean = getMean(sample, dim, iweight)
502 8142 : cor_ref = getCorRef(mean, sample, dim, real(iweight, TKC))
503 :
504 6360 : cor = getFilled(ZERO, ndim, ndim)
505 400 : call setCor(cor, subsetr, mean, sample, dim, iweight, iweisum)
506 400 : call report(__LINE__, subsetr, SK_"integer-weighted")
507 :
508 6360 : cor = getFilled(ZERO, ndim, ndim)
509 1595 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim, iweight, iweisum)
510 400 : call report(__LINE__, subsetr, SK_"integer-weighted shifted")
511 :
512 : ! real weighted
513 :
514 1995 : mean = getMean(sample, dim, rweight)
515 6360 : cor_ref = getCorRef(mean, sample, dim, rweight)
516 :
517 6360 : cor = getFilled(ZERO, ndim, ndim)
518 400 : call setCor(cor, subsetr, mean, sample, dim, rweight, rweisum)
519 400 : call report(__LINE__, subsetr, SK_"real-weighted")
520 :
521 6360 : cor = getFilled(ZERO, ndim, ndim)
522 1595 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim, rweight, rweisum)
523 400 : call report(__LINE__, subsetr, SK_"real-weighted shifted")
524 :
525 : ! unweighted
526 :
527 1995 : mean = getMean(sample, dim)
528 6360 : cor_ref = getCorRef(mean, sample, dim)
529 :
530 6360 : cor = getFilled(ZERO, ndim, ndim)
531 400 : call setCor(cor, subsetr, mean, sample, dim)
532 400 : call report(__LINE__, subsetr, SK_"unweighted")
533 :
534 6360 : cor = getFilled(ZERO, ndim, ndim)
535 1595 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim)
536 400 : call report(__LINE__, subsetr, SK_"unweighted shifted")
537 :
538 : end block
539 :
540 : ! test sample uppDia interface.
541 :
542 : block
543 :
544 : type(uppDia_type), parameter :: subsetr = uppDia_type()
545 400 : dim = getChoice([1, 2])
546 400 : ndim = getUnifRand(1_IK, 5_IK)
547 400 : if (dim == 2) then
548 3946 : sample = getUnifRand(ONES, TWOS, ndim, nsam)
549 : else
550 3859 : sample = getUnifRand(ONES, TWOS, nsam, ndim)
551 : end if
552 2582 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
553 2582 : iweight = getUnifRand(1_IK, 9_IK, nsam)
554 2182 : iweisum = sum(iweight)
555 2182 : rweisum = sum(rweight)
556 :
557 : ! integer weighted
558 :
559 2022 : mean = getMean(sample, dim, iweight)
560 8310 : cor_ref = getCorRef(mean, sample, dim, real(iweight, TKC))
561 :
562 6528 : cor = getFilled(ZERO, ndim, ndim)
563 400 : call setCor(cor, subsetr, mean, sample, dim, iweight, iweisum)
564 400 : call report(__LINE__, subsetr, SK_"integer-weighted")
565 :
566 6528 : cor = getFilled(ZERO, ndim, ndim)
567 1622 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim, iweight, iweisum)
568 400 : call report(__LINE__, subsetr, SK_"integer-weighted shifted")
569 :
570 : ! real weighted
571 :
572 2022 : mean = getMean(sample, dim, rweight)
573 6528 : cor_ref = getCorRef(mean, sample, dim, rweight)
574 :
575 6528 : cor = getFilled(ZERO, ndim, ndim)
576 400 : call setCor(cor, subsetr, mean, sample, dim, rweight, rweisum)
577 400 : call report(__LINE__, subsetr, SK_"real-weighted")
578 :
579 6528 : cor = getFilled(ZERO, ndim, ndim)
580 1622 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim, rweight, rweisum)
581 400 : call report(__LINE__, subsetr, SK_"real-weighted shifted")
582 :
583 : ! unweighted
584 :
585 2022 : mean = getMean(sample, dim)
586 6528 : cor_ref = getCorRef(mean, sample, dim)
587 :
588 6528 : cor = getFilled(ZERO, ndim, ndim)
589 400 : call setCor(cor, subsetr, mean, sample, dim)
590 400 : call report(__LINE__, subsetr, SK_"unweighted")
591 :
592 6528 : cor = getFilled(ZERO, ndim, ndim)
593 1622 : call setCor(cor, subsetr, getShifted(sample, dim, -mean), dim)
594 400 : call report(__LINE__, subsetr, SK_"unweighted shifted")
595 :
596 : end block
597 :
598 : ! test sample XY interface.
599 :
600 8 : block
601 :
602 : type(lowDia_type), parameter :: subsetr = lowDia_type()
603 400 : dim = 1_IK
604 400 : ndim = 2_IK
605 1600 : mean = getFilled(ZERO, ndim)
606 3200 : cor_ref = getFilled(ZERO, ndim, ndim)
607 5164 : sample = getUnifRand(ONES, TWOS, nsam, ndim)
608 2582 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
609 2582 : iweight = getUnifRand(1_IK, 9_IK, nsam)
610 2182 : iweisum = sum(iweight)
611 2182 : rweisum = sum(rweight)
612 :
613 : ! integer weighted
614 :
615 1600 : mean = getMean(sample, dim, iweight)
616 4982 : cor_ref = getCorRef(mean, sample, dim, real(iweight, TKC))
617 :
618 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
619 400 : call setCor(cor(1,2), mean, sample(:,1), sample(:,2), iweight, iweisum)
620 400 : call report(__LINE__, uppLowDia, SK_"integer-weighted")
621 :
622 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
623 3964 : call setCor(cor(1,2), sample(:,1) - mean(1), sample(:,2) - mean(2), iweight, iweisum)
624 400 : call report(__LINE__, uppLowDia, SK_"integer-weighted")
625 :
626 : ! real weighted
627 :
628 1600 : mean = getMean(sample, dim, rweight)
629 3200 : cor_ref = getCorRef(mean, sample, dim, rweight)
630 :
631 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
632 400 : call setCor(cor(1,2), mean, sample(:,1), sample(:,2), rweight, rweisum)
633 400 : call report(__LINE__, uppLowDia, SK_"real-weighted")
634 :
635 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
636 3964 : call setCor(cor(1,2), sample(:,1) - mean(1), sample(:,2) - mean(2), rweight, rweisum)
637 400 : call report(__LINE__, uppLowDia, SK_"real-weighted")
638 :
639 : ! unweighted
640 :
641 1600 : mean = getMean(sample, dim)
642 3200 : cor_ref = getCorRef(mean, sample, dim)
643 :
644 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
645 400 : call setCor(cor(1,2), mean, sample(:,1), sample(:,2))
646 400 : call report(__LINE__, uppLowDia, SK_"unweighted")
647 :
648 6400 : cor = getMatInit([ndim, ndim], uppLowDia, ZERO, ZERO, ONE)
649 3964 : call setCor(cor(1,2), sample(:,1) - mean(1), sample(:,2) - mean(2))
650 400 : call report(__LINE__, uppLowDia, SK_"unweighted")
651 :
652 : end block
653 :
654 : end do
655 :
656 : contains
657 :
658 6400 : subroutine reportCFC(line, stdinv)
659 : real(TKC), intent(in), optional :: stdinv(:)
660 : integer, intent(in) :: line
661 147656 : cdiff = abs(cor - cor_ref)
662 95568 : assertion = assertion .and. all(cdiff < ctol)
663 6400 : if (test%traceable .and. .not. assertion) then
664 : ! LCOV_EXCL_START
665 : call test%disp%skip()
666 : call test%disp%show("ndim")
667 : call test%disp%show( ndim )
668 : call test%disp%show("present(stdinv)")
669 : call test%disp%show( present(stdinv) )
670 : if (present(stdinv)) then
671 : call test%disp%show("stdinv")
672 : call test%disp%show( stdinv )
673 : end if
674 : call test%disp%show("cov")
675 : call test%disp%show( cov )
676 : call test%disp%show("cor")
677 : call test%disp%show( cor )
678 : call test%disp%show("cor_ref")
679 : call test%disp%show( cor_ref )
680 : call test%disp%show("cdiff")
681 : call test%disp%show( cdiff )
682 : call test%disp%skip()
683 : ! LCOV_EXCL_STOP
684 : end if
685 6400 : call test%assert(assertion, SK_"The `cor` must be correctly computed from the specified `cor` and `stdinv`.", int(line, IK))
686 6400 : end subroutine
687 :
688 3600 : PURE function getCorRef(mean, sample, dim, weight) result(cor)
689 : ! returns the full correlation matrix of the sample.
690 : integer(IK), intent(in) :: dim
691 : real(TKC), intent(in), optional :: weight(:)
692 : TYPE_OF_SAMPLE, intent(in), contiguous :: sample(:,:), mean(:)
693 : TYPE_OF_SAMPLE :: cor(size(sample, 3 - dim, IK), size(sample, 3 - dim, IK))
694 7200 : TYPE_OF_SAMPLE :: sampleShifted(size(sample, 1, IK), size(sample, 2, IK))
695 : integer(IK) :: idim, jdim, ndim
696 7200 : real(TKC) :: stdinv(size(sample, 3 - dim, IK))
697 3600 : ndim = size(sample, 3 - dim, IK)
698 13251 : sampleShifted = getShifted(sample, dim, -mean)
699 3600 : if (present(weight)) then
700 23784 : call setCov(cor, uppDia, sampleShifted, dim, weight, sum(weight))
701 : else
702 1200 : call setCov(cor, uppDia, sampleShifted, dim)
703 : end if
704 13251 : do jdim = 1, ndim
705 9651 : stdinv(jdim) = 1._TKC / sqrt(real(cor(jdim, jdim), TKC))
706 9651 : cor(jdim, jdim) = 1._TKC
707 24132 : do idim = jdim - 1, 1, -1
708 10881 : cor(idim, jdim) = cor(idim, jdim) * stdinv(idim) * stdinv(jdim)
709 20532 : cor(jdim, idim) = GET_CONJG(cor(idim, jdim))
710 : end do
711 : end do
712 3600 : end function
713 :
714 7200 : subroutine report(line, subsetr, this)
715 : integer, intent(in) :: line
716 : class(*), intent(in) :: subsetr
717 : character(*, SK), intent(in) :: this
718 7200 : if (same_type_as(subsetr, uppDia)) then
719 2400 : call setMatCopy(cor, rdpack, cor, rdpack, uppDia, transHerm)
720 4800 : elseif (same_type_as(subsetr, lowDia)) then
721 2400 : call setMatCopy(cor, rdpack, cor, rdpack, lowDia, transHerm)
722 2400 : elseif (same_type_as(subsetr, uppLowDia)) then
723 2400 : cor(2,1) = GET_CONJG(cor(1,2))
724 : else
725 0 : error stop "Unrecognized subset."
726 : end if
727 138120 : cdiff = abs(cor - cor_ref)
728 89328 : assertion = assertion .and. all(cdiff < ctol)
729 7200 : if (test%traceable .and. .not. assertion) then
730 : ! LCOV_EXCL_START
731 : call test%disp%skip()
732 : call test%disp%show("[ndim, nsam, dim, shape(sample, IK)]")
733 : call test%disp%show( [ndim, nsam, dim, shape(sample, IK)] )
734 : call test%disp%show("sample")
735 : call test%disp%show( sample )
736 : call test%disp%show("iweight")
737 : call test%disp%show( iweight )
738 : call test%disp%show("rweight")
739 : call test%disp%show( rweight )
740 : call test%disp%show("cor_ref")
741 : call test%disp%show( cor_ref )
742 : call test%disp%show("cor")
743 : call test%disp%show( cor )
744 : call test%disp%show("cdiff")
745 : call test%disp%show( cdiff )
746 : call test%disp%skip()
747 : ! LCOV_EXCL_STOP
748 : end if
749 7200 : call test%assert(assertion, SK_"The `cor` must be computed correctly for "//this//SK_" sample.", int(line, IK))
750 7200 : end subroutine
751 :
752 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
753 : #elif getRho_ENABLED && (XY_ENABLED || D0_ENABLED)
754 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
755 :
756 : integer(IK) :: itry, nsam
757 : real(RK), allocatable :: rweight(:)
758 : integer(IK), allocatable :: iweight(:)
759 : real(RK), allocatable :: frankx(:), franky(:)
760 : #if D0_ENABLED
761 1 : TYPE_OF_SAMPLE, allocatable :: x, y
762 : #elif XY_ENABLED
763 1 : TYPE_OF_SAMPLE, allocatable :: x(:), y(:)
764 : #else
765 : #error "Unrecognized interface."
766 : #endif
767 : real(RK) :: rho, rho_ref, rdiff
768 12 : assertion = .true._LK
769 666 : do itry = 1, 50
770 :
771 600 : nsam = getUnifRand(20_IK, 30_IK)
772 : #if D0_ENABLED && SK_ENABLED
773 2552 : x = getUnifRand(repeat(lb, nsam), repeat(ub, nsam))
774 2552 : y = getUnifRand(repeat(lb, nsam), repeat(ub, nsam))
775 : #elif XY_ENABLED && (BSSK_ENABLED || PSSK_ENABLED)
776 2498 : if (allocated(x)) deallocate(x); allocate(x(nsam))
777 2498 : if (allocated(y)) deallocate(y); allocate(y(nsam))
778 : block
779 : integer(IK) :: isam
780 1287 : do isam = 1, nsam
781 1237 : x(isam)%val = getUnifRand(lb, ub)
782 1287 : y(isam)%val = getUnifRand(lb, ub)
783 : end do
784 : end block
785 : #elif XY_ENABLED
786 13504 : x = getUnifRand(lb, ub, nsam)
787 13504 : y = getUnifRand(lb, ub, nsam)
788 : #else
789 : #error "Unrecognized interface."
790 : #endif
791 16192 : iweight = getUnifRand(1_IK, 9_IK, nsam)
792 16192 : rweight = getUnifRand(1._RK, 9._RK, nsam)
793 :
794 : ! integer weighted
795 :
796 15592 : rho_ref = getRhoCoef(real(iweight, RK))
797 :
798 600 : rho = getRho(x, y, iweight)
799 600 : call report(__LINE__, SK_"integer-weighted")
800 :
801 : ! real weighted
802 :
803 600 : rho_ref = getRhoCoef(rweight)
804 :
805 600 : rho = getRho(x, y, rweight)
806 600 : call report(__LINE__, SK_"real-weighted")
807 :
808 : ! unweighted
809 :
810 600 : rho_ref = getRhoCoef()
811 :
812 600 : rho = getRho(x, y)
813 612 : call report(__LINE__, SK_"unweighted")
814 :
815 : end do
816 :
817 : contains
818 :
819 1800 : function getRhoCoef(weight) result(rho)
820 : ! returns the full correlation matrix of the sample.
821 : real(RK), intent(in), optional :: weight(:)
822 : real(RK) :: rho
823 48576 : frankx = getRankFractional(x)
824 48576 : franky = getRankFractional(y)
825 1800 : if (present(weight)) then
826 31184 : rho = getCor(frankx, franky, weight)
827 : else
828 600 : rho = getCor(frankx, franky)
829 : end if
830 1800 : end function
831 :
832 1800 : subroutine report(line, this)
833 : integer, intent(in) :: line
834 : character(*, SK), intent(in) :: this
835 1800 : rdiff = abs(rho - rho_ref)
836 1800 : assertion = assertion .and. rdiff < rtol
837 1800 : if (test%traceable .and. .not. assertion) then
838 : ! LCOV_EXCL_START
839 : call test%disp%skip()
840 : call test%disp%show("nsam")
841 : call test%disp%show( nsam )
842 : call test%disp%show("x")
843 : call test%disp%show( x , deliml = SK_"""" )
844 : call test%disp%show("y" )
845 : call test%disp%show( y , deliml = SK_"""" )
846 : call test%disp%show("iweight")
847 : call test%disp%show( iweight )
848 : call test%disp%show("rweight")
849 : call test%disp%show( rweight )
850 : call test%disp%show("rho_ref")
851 : call test%disp%show( rho_ref )
852 : call test%disp%show("rho")
853 : call test%disp%show( rho )
854 : call test%disp%show("rdiff")
855 : call test%disp%show( rdiff )
856 : call test%disp%skip()
857 : ! LCOV_EXCL_STOP
858 : end if
859 1800 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
860 1800 : end subroutine
861 :
862 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
863 : #elif getRho_ENABLED && D2_ENABLED
864 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
865 :
866 : real(RK), allocatable :: rweight(:)
867 : integer(IK), allocatable :: iweight(:)
868 : integer(IK) :: itry, nsam
869 : real(RK), allocatable :: rho(:,:), rho_ref(:,:), rdiff(:,:)
870 1 : TYPE_OF_SAMPLE, allocatable :: sample(:,:)
871 : integer(IK) :: ndim, dim
872 11 : assertion = .true._LK
873 :
874 584 : do itry = 1, 50
875 :
876 : ! Test Rho: The strategy is to construct a sample, compute its correlation and compare it against the reference computed here.
877 :
878 550 : nsam = getUnifRand(20_IK, 30_IK)
879 :
880 : ! test sample lowDia interface.
881 :
882 11 : block
883 :
884 550 : dim = getChoice([1, 2])
885 550 : ndim = getUnifRand(1_IK, 5_IK)
886 550 : if (dim == 2) then
887 : #if PSSK_ENABLED
888 : sample = css_pdt(getUnifRand(lb, ub, ndim, nsam))
889 : #elif BSSK_ENABLED
890 8035 : sample = css_type(getUnifRand(lb, ub, ndim, nsam))
891 : #else
892 25930 : sample = getUnifRand(lb, ub, ndim, nsam)
893 : #endif
894 : else
895 : #if PSSK_ENABLED
896 : sample = css_pdt(getUnifRand(lb, ub, nsam, ndim))
897 : #elif BSSK_ENABLED
898 7769 : sample = css_type(getUnifRand(lb, ub, nsam, ndim))
899 : #else
900 21573 : sample = getUnifRand(lb, ub, nsam, ndim)
901 : #endif
902 : end if
903 14945 : iweight = getUnifRand(1_IK, 9_IK, nsam)
904 14945 : rweight = getUnifRand(1., 9., nsam)
905 :
906 : ! integer weighted
907 :
908 23075 : rho_ref = getRhoRef(real(iweight, RK))
909 9230 : rho = getRho(sample, dim, iweight)
910 550 : call report(__LINE__, SK_"integer-weighted")
911 :
912 : ! real weighted
913 :
914 9230 : rho_ref = getRhoRef(rweight)
915 9230 : rho = getRho(sample, dim, rweight)
916 550 : call report(__LINE__, SK_"real-weighted")
917 :
918 : ! unweighted
919 :
920 9230 : rho_ref = getRhoRef()
921 9230 : rho = getRho(sample, dim)
922 550 : call report(__LINE__, SK_"unweighted")
923 :
924 : end block
925 :
926 : end do
927 :
928 : contains
929 :
930 1650 : function getRhoRef(weight) result(rho)
931 : ! returns the full correlation matrix of the sample.
932 : real(RK), intent(in), optional :: weight(:)
933 : real(RK) :: rho(size(sample, 3 - dim, IK), size(sample, 3 - dim, IK))
934 1650 : real(RK) :: frank(size(sample, 1, IK), size(sample, 2, IK))
935 : integer(IK) :: idim, ndim, nsam
936 : ndim = size(sample, 3 - dim, IK)
937 : nsam = size(sample, dim, IK)
938 1650 : if (dim == 2_IK) then
939 3327 : do idim = 1, ndim
940 129507 : frank(idim, :) = getRankFractional(sample(idim, :))
941 : end do
942 : else
943 3462 : do idim = 1, ndim
944 3462 : frank(:, idim) = getRankFractional(sample(:, idim))
945 : end do
946 : end if
947 1650 : if (present(weight)) then
948 28790 : rho = getCor(frank, dim, weight)
949 : else
950 550 : rho = getCor(frank, dim)
951 : end if
952 1650 : end function
953 :
954 1650 : subroutine report(line, this)
955 : integer, intent(in) :: line
956 : character(*, SK), intent(in) :: this
957 27690 : rdiff = abs(rho - rho_ref)
958 26040 : assertion = assertion .and. all(rdiff < rtol)
959 1650 : if (test%traceable .and. .not. assertion) then
960 : ! LCOV_EXCL_START
961 : call test%disp%skip()
962 : call test%disp%show("[ndim, nsam, dim, GET_SHAPE(sample, IK)]")
963 : call test%disp%show( [ndim, nsam, dim, GET_SHAPE(sample, IK)] )
964 : call test%disp%show("sample")
965 : call test%disp%show( sample )
966 : call test%disp%show("iweight")
967 : call test%disp%show( iweight )
968 : call test%disp%show("rweight")
969 : call test%disp%show( rweight )
970 : call test%disp%show("rho_ref")
971 : call test%disp%show( rho_ref )
972 : call test%disp%show("rho")
973 : call test%disp%show( rho )
974 : call test%disp%show("rdiff")
975 : call test%disp%show( rdiff )
976 : call test%disp%skip()
977 : ! LCOV_EXCL_STOP
978 : end if
979 1650 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
980 1650 : end subroutine
981 :
982 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
983 : #elif setRho_ENABLED && (XY_ENABLED || D0_ENABLED)
984 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
985 :
986 : integer(IK) :: itry, nsam
987 : real(RK), allocatable :: rweight(:)
988 : integer(IK), allocatable :: iweight(:)
989 : real(RK), allocatable :: frankx(:), franky(:), frankx_ref(:), franky_ref(:)
990 : #if D0_ENABLED
991 1 : TYPE_OF_SAMPLE, allocatable :: x, y
992 : #elif XY_ENABLED
993 1 : TYPE_OF_SAMPLE, allocatable :: x(:), y(:)
994 : #else
995 : #error "Unrecognized interface."
996 : #endif
997 : real(RK) :: rho, rho_ref, rdiff
998 12 : assertion = .true._LK
999 668 : do itry = 1, 50
1000 :
1001 600 : nsam = getUnifRand(20_IK, 30_IK)
1002 : #if D0_ENABLED && SK_ENABLED
1003 2536 : x = getUnifRand(repeat(lb, nsam), repeat(ub, nsam))
1004 2536 : y = getUnifRand(repeat(lb, nsam), repeat(ub, nsam))
1005 : #elif XY_ENABLED && (BSSK_ENABLED || PSSK_ENABLED)
1006 2481 : if (allocated(x)) deallocate(x); allocate(x(nsam))
1007 2481 : if (allocated(y)) deallocate(y); allocate(y(nsam))
1008 : block
1009 : integer(IK) :: isam
1010 1279 : do isam = 1, nsam
1011 1229 : x(isam)%val = getUnifRand(lb, ub)
1012 1279 : y(isam)%val = getUnifRand(lb, ub)
1013 : end do
1014 : end block
1015 : #elif XY_ENABLED
1016 13469 : x = getUnifRand(lb, ub, nsam)
1017 13469 : y = getUnifRand(lb, ub, nsam)
1018 : #else
1019 : #error "Unrecognized interface."
1020 : #endif
1021 16141 : frankx = getFilled(0._RK, nsam)
1022 16141 : franky = getFilled(0._RK, nsam)
1023 16141 : iweight = getUnifRand(1_IK, 9_IK, nsam)
1024 16141 : rweight = getUnifRand(1._RK, 9._RK, nsam)
1025 :
1026 : ! integer weighted
1027 :
1028 15541 : rho_ref = getRhoCoef(real(iweight, RK))
1029 :
1030 : call setRho(rho, frankx, franky, x, y, iweight)
1031 600 : call report(__LINE__, SK_"integer-weighted")
1032 :
1033 : ! real weighted
1034 :
1035 600 : rho_ref = getRhoCoef(rweight)
1036 :
1037 : call setRho(rho, frankx, franky, x, y, rweight)
1038 600 : call report(__LINE__, SK_"real-weighted")
1039 :
1040 : ! unweighted
1041 :
1042 600 : rho_ref = getRhoCoef()
1043 :
1044 : call setRho(rho, frankx, franky, x, y)
1045 1812 : call report(__LINE__, SK_"unweighted")
1046 :
1047 : end do
1048 :
1049 : contains
1050 :
1051 1800 : function getRhoCoef(weight) result(rho)
1052 : ! returns the full correlation matrix of the sample.
1053 : real(RK), intent(in), optional :: weight(:)
1054 : real(RK) :: rho
1055 48423 : frankx_ref = getRankFractional(x)
1056 48423 : franky_ref = getRankFractional(y)
1057 1800 : if (present(weight)) then
1058 31082 : rho = getCor(frankx_ref, franky_ref, weight)
1059 : else
1060 600 : rho = getCor(frankx_ref, franky_ref)
1061 : end if
1062 1800 : end function
1063 :
1064 1800 : subroutine report(line, this)
1065 : integer, intent(in) :: line
1066 : character(*, SK), intent(in) :: this
1067 1800 : rdiff = abs(rho - rho_ref)
1068 1800 : assertion = assertion .and. rdiff < rtol
1069 46623 : assertion = assertion .and. all(abs(frankx - frankx_ref) < rtol)
1070 46623 : assertion = assertion .and. all(abs(franky - franky_ref) < rtol)
1071 1800 : if (test%traceable .and. .not. assertion) then
1072 : ! LCOV_EXCL_START
1073 : call test%disp%skip()
1074 : call test%disp%show("nsam")
1075 : call test%disp%show( nsam )
1076 : call test%disp%show("x")
1077 : call test%disp%show( x , deliml = SK_"""" )
1078 : call test%disp%show("y" )
1079 : call test%disp%show( y , deliml = SK_"""" )
1080 : call test%disp%show("iweight")
1081 : call test%disp%show( iweight )
1082 : call test%disp%show("rweight")
1083 : call test%disp%show( rweight )
1084 : call test%disp%show("frankx")
1085 : call test%disp%show( frankx )
1086 : call test%disp%show("frankx_ref")
1087 : call test%disp%show( frankx_ref )
1088 : call test%disp%show("franky")
1089 : call test%disp%show( franky )
1090 : call test%disp%show("franky_ref")
1091 : call test%disp%show( franky_ref )
1092 : call test%disp%show("rho_ref")
1093 : call test%disp%show( rho_ref )
1094 : call test%disp%show("rho")
1095 : call test%disp%show( rho )
1096 : call test%disp%show("rdiff")
1097 : call test%disp%show( rdiff )
1098 : call test%disp%skip()
1099 : ! LCOV_EXCL_STOP
1100 : end if
1101 : ! Three tests, so three calls.
1102 1800 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
1103 1800 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
1104 1800 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
1105 1800 : end subroutine
1106 :
1107 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
1108 : #elif setRho_ENABLED && D2_ENABLED
1109 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
1110 :
1111 : real(RK), allocatable :: rweight(:)
1112 : integer(IK), allocatable :: iweight(:)
1113 : integer(IK) :: itry, nsam
1114 : real(RK), allocatable :: frank(:,:), frank_ref(:,:), rho(:,:), rho_ref(:,:), rdiff(:,:)
1115 1 : TYPE_OF_SAMPLE, allocatable :: sample(:,:)
1116 : integer(IK) :: ndim, dim
1117 11 : assertion = .true._LK
1118 :
1119 592 : do itry = 1, 50
1120 :
1121 : ! Test Rho: The strategy is to construct a sample, compute its correlation and compare it against the reference computed here.
1122 :
1123 550 : nsam = getUnifRand(20_IK, 30_IK)
1124 :
1125 : ! test sample lowDia interface.
1126 :
1127 : block
1128 :
1129 : type(lowDia_type), parameter :: subsetr = lowDia_type()
1130 550 : dim = getChoice([1, 2])
1131 550 : ndim = getUnifRand(1_IK, 5_IK)
1132 550 : if (dim == 2) then
1133 : #if PSSK_ENABLED
1134 : sample = css_pdt(getUnifRand(lb, ub, ndim, nsam))
1135 : #elif BSSK_ENABLED
1136 7993 : sample = css_type(getUnifRand(lb, ub, ndim, nsam))
1137 : #else
1138 25868 : sample = getUnifRand(lb, ub, ndim, nsam)
1139 : #endif
1140 : else
1141 : #if PSSK_ENABLED
1142 : sample = css_pdt(getUnifRand(lb, ub, nsam, ndim))
1143 : #elif BSSK_ENABLED
1144 6074 : sample = css_type(getUnifRand(lb, ub, nsam, ndim))
1145 : #else
1146 19391 : sample = getUnifRand(lb, ub, nsam, ndim)
1147 : #endif
1148 : end if
1149 49252 : frank = getFilled(0._RK, size(sample, 1, IK), size(sample, 2, IK))
1150 50336 : frank_ref = frank
1151 14779 : iweight = getUnifRand(1_IK, 9_IK, nsam)
1152 14779 : rweight = getUnifRand(1., 9., nsam)
1153 :
1154 : ! integer weighted
1155 :
1156 22389 : rho_ref = getRhoRef(real(iweight, RK))
1157 :
1158 8710 : rho = getFilled(0._RK, ndim, ndim)
1159 550 : call setRho(rho, subsetr, frank, sample, dim, iweight)
1160 550 : call report(__LINE__, subsetr, SK_"integer-weighted")
1161 :
1162 : ! real weighted
1163 :
1164 8710 : rho_ref = getRhoRef(rweight)
1165 :
1166 8710 : rho = getFilled(0._RK, ndim, ndim)
1167 550 : call setRho(rho, subsetr, frank, sample, dim, rweight)
1168 550 : call report(__LINE__, subsetr, SK_"real-weighted")
1169 :
1170 : ! unweighted
1171 :
1172 8710 : rho_ref = getRhoRef()
1173 :
1174 8710 : rho = getFilled(0._RK, ndim, ndim)
1175 550 : call setRho(rho, subsetr, frank, sample, dim)
1176 550 : call report(__LINE__, subsetr, SK_"unweighted")
1177 :
1178 : end block
1179 :
1180 : ! test sample uppDia interface.
1181 :
1182 11 : block
1183 :
1184 : type(uppDia_type), parameter :: subsetr = uppDia_type()
1185 550 : dim = getChoice([1, 2])
1186 550 : ndim = getUnifRand(1_IK, 5_IK)
1187 550 : if (dim == 2) then
1188 : #if PSSK_ENABLED
1189 : sample = css_pdt(getUnifRand(lb, ub, ndim, nsam))
1190 : #elif BSSK_ENABLED
1191 7174 : sample = css_type(getUnifRand(lb, ub, ndim, nsam))
1192 : #else
1193 25801 : sample = getUnifRand(lb, ub, ndim, nsam)
1194 : #endif
1195 : else
1196 : #if PSSK_ENABLED
1197 : sample = css_pdt(getUnifRand(lb, ub, nsam, ndim))
1198 : #elif BSSK_ENABLED
1199 6085 : sample = css_type(getUnifRand(lb, ub, nsam, ndim))
1200 : #else
1201 20168 : sample = getUnifRand(lb, ub, nsam, ndim)
1202 : #endif
1203 : end if
1204 49926 : frank = getFilled(0._RK, size(sample, 1, IK), size(sample, 2, IK))
1205 50940 : frank_ref = frank
1206 14779 : iweight = getUnifRand(1_IK, 9_IK, nsam)
1207 14779 : rweight = getUnifRand(1., 9., nsam)
1208 :
1209 : ! integer weighted
1210 :
1211 22569 : rho_ref = getRhoRef(real(iweight, RK))
1212 :
1213 8890 : rho = getFilled(0._RK, ndim, ndim)
1214 550 : call setRho(rho, subsetr, frank, sample, dim, iweight)
1215 550 : call report(__LINE__, subsetr, SK_"integer-weighted")
1216 :
1217 : ! real weighted
1218 :
1219 8890 : rho_ref = getRhoRef(rweight)
1220 :
1221 8890 : rho = getFilled(0._RK, ndim, ndim)
1222 550 : call setRho(rho, subsetr, frank, sample, dim, rweight)
1223 550 : call report(__LINE__, subsetr, SK_"real-weighted")
1224 :
1225 : ! unweighted
1226 :
1227 8890 : rho_ref = getRhoRef()
1228 :
1229 8890 : rho = getFilled(0._RK, ndim, ndim)
1230 550 : call setRho(rho, subsetr, frank, sample, dim)
1231 550 : call report(__LINE__, subsetr, SK_"unweighted")
1232 :
1233 : end block
1234 :
1235 : end do
1236 :
1237 : contains
1238 :
1239 3300 : function getRhoRef(weight) result(rho)
1240 : ! returns the full correlation matrix of the sample.
1241 : real(RK), intent(in), optional :: weight(:)
1242 : real(RK) :: rho(size(sample, 3 - dim, IK), size(sample, 3 - dim, IK))
1243 : integer(IK) :: idim, ndim, nsam
1244 : ndim = size(sample, 3 - dim, IK)
1245 : nsam = size(sample, dim, IK)
1246 3300 : if (dim == 2_IK) then
1247 6660 : do idim = 1, ndim
1248 254634 : frank_ref(idim, :) = getRankFractional(sample(idim, :))
1249 : end do
1250 : else
1251 6486 : do idim = 1, ndim
1252 127290 : frank_ref(:, idim) = getRankFractional(sample(:, idim))
1253 : end do
1254 : end if
1255 3300 : if (present(weight)) then
1256 56916 : rho = getCor(frank_ref, dim, weight)
1257 : else
1258 1100 : rho = getCor(frank_ref, dim)
1259 : end if
1260 3300 : end function
1261 :
1262 3300 : subroutine report(line, subsetr, this)
1263 : integer, intent(in) :: line
1264 : class(*), intent(in) :: subsetr
1265 : character(*, SK), intent(in) :: this
1266 3300 : if (same_type_as(subsetr, uppDia)) then
1267 1650 : call setMatCopy(rho, rdpack, rho, rdpack, uppDia, transHerm)
1268 1650 : elseif (same_type_as(subsetr, lowDia)) then
1269 1650 : call setMatCopy(rho, rdpack, rho, rdpack, lowDia, transHerm)
1270 0 : elseif (same_type_as(subsetr, uppLowDia)) then
1271 0 : rho(2,1) = rho(1,2)
1272 : else
1273 0 : error stop "Unrecognized subset."
1274 : end if
1275 52800 : rdiff = abs(rho - rho_ref)
1276 49500 : assertion = assertion .and. all(rdiff < rtol)
1277 294234 : assertion = assertion .and. all(abs(frank - frank_ref) < rtol)
1278 3300 : if (test%traceable .and. .not. assertion) then
1279 : ! LCOV_EXCL_START
1280 : call test%disp%skip()
1281 : call test%disp%show("[ndim, nsam, dim, GET_SHAPE(sample, IK)]")
1282 : call test%disp%show( [ndim, nsam, dim, GET_SHAPE(sample, IK)] )
1283 : call test%disp%show("sample")
1284 : call test%disp%show( sample )
1285 : call test%disp%show("iweight")
1286 : call test%disp%show( iweight )
1287 : call test%disp%show("rweight")
1288 : call test%disp%show( rweight )
1289 : call test%disp%show("frank_ref")
1290 : call test%disp%show( frank_ref )
1291 : call test%disp%show("frank")
1292 : call test%disp%show( frank )
1293 : call test%disp%show("rho_ref")
1294 : call test%disp%show( rho_ref )
1295 : call test%disp%show("rho")
1296 : call test%disp%show( rho )
1297 : call test%disp%show("rdiff")
1298 : call test%disp%show( rdiff )
1299 : call test%disp%skip()
1300 : ! LCOV_EXCL_STOP
1301 : end if
1302 : ! two tests, two reports.
1303 3300 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
1304 3300 : call test%assert(assertion, SK_"The `rho` and fractional ranks must be computed correctly for "//this//SK_" sample.", int(line, IK))
1305 3300 : end subroutine
1306 :
1307 : !%%%%%%%%%%%%%%%%%%
1308 : #elif setCordance_ENABLED
1309 : !%%%%%%%%%%%%%%%%%%
1310 :
1311 : assertion = .true.
1312 :
1313 : #else
1314 : !%%%%%%%%%%%%%%%%%%%%%%%%
1315 : #error "Unrecognized interface."
1316 : !%%%%%%%%%%%%%%%%%%%%%%%%
1317 : #endif
1318 : #undef TYPE_OF_SAMPLE
1319 : #undef GET_CONJG
1320 : #undef GET_SHAPE
|
