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_sampleCov](@ref pm_sampleCov).
19 : !>
20 : !> \fintest
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 5:03 PM, August 11, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : real(TKC), parameter :: rtol = epsilon(1._TKC) * 100
28 : ! Define the sample type.
29 : #if CK_ENABLED
30 : #define GET_CONJG(X)conjg(X)
31 : #define TYPE_OF_SAMPLE complex(TKC)
32 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 1._TKC), TWO = 2 * (1._TKC, 1._TKC), ctol = (rtol, rtol), RONE = (1._TKC, 0._TKC)
33 : #elif RK_ENABLED
34 : #define GET_CONJG(X)X
35 : #define TYPE_OF_SAMPLE real(TKC)
36 : real(TKC), parameter :: ZERO = 0._TKC, ONE = 1._TKC, RONE = 1._TKC, TWO = 2._TKC, ctol = rtol
37 : #else
38 : #error "Unrecognized interface."
39 : #endif
40 : !%%%%%%%%%%%%%
41 : #if getCov_ENABLED
42 : !%%%%%%%%%%%%%
43 :
44 : real(TKC) :: rweisqs
45 : real(TKC) :: rweisum
46 : integer(IK) :: iweisum
47 : real(TKC), allocatable :: rweight(:)
48 : integer(IK), allocatable :: iweight(:)
49 : integer(IK) :: itry, nsam, ndim, dim
50 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:), meang(:)
51 : TYPE_OF_SAMPLE, allocatable :: cov(:,:), cor(:,:), cov_ref(:,:), cdiff(:,:)
52 : real(TKC), allocatable :: std(:)
53 8 : assertion = .true._LK
54 :
55 408 : do itry = 1, 50
56 :
57 : ! test CFC subsetr = lowDia.
58 :
59 : block
60 :
61 : type(lowDia_type), parameter :: subsetr = lowDia_type()
62 400 : ndim = getUnifRand(1_IK, 5_IK)
63 6548 : cov_ref = getFilled(ZERO, ndim, ndim)
64 2022 : std = getUnifRand(1._TKC, 2._TKC, ndim)
65 6548 : cor = getUnifRand(-ONE, ONE, ndim, ndim)
66 400 : call setMatInit(cor, getSubSymm(subsetr), ZERO, RONE)
67 400 : call setCov(cov_ref, subsetr, cor, subsetr, std)
68 400 : call setMatCopy(cov_ref, rdpack, cov_ref, rdpack, subsetr, transHerm)
69 :
70 12296 : cov = getCov(cor, subsetr, std)
71 400 : call reportCFC(__LINE__)
72 :
73 : end block
74 :
75 : ! test CFC subsetr = uppDia.
76 :
77 : block
78 :
79 : type(uppDia_type), parameter :: subsetr = uppDia_type()
80 400 : ndim = getUnifRand(1_IK, 5_IK)
81 6496 : cov_ref = getFilled(ZERO, ndim, ndim)
82 2020 : std = getUnifRand(1._TKC, 2._TKC, ndim)
83 6496 : cor = getUnifRand(-ONE, ONE, ndim, ndim)
84 400 : call setMatInit(cor, getSubSymm(subsetr), ZERO, RONE)
85 400 : call setCov(cov_ref, subsetr, cor, subsetr, std)
86 400 : call setMatCopy(cov_ref, rdpack, cov_ref, rdpack, subsetr, transHerm)
87 :
88 12192 : cov = getCov(cor, subsetr, std)
89 400 : call reportCFC(__LINE__)
90 :
91 : end block
92 :
93 400 : nsam = getUnifRand(2_IK, 7_IK)
94 :
95 : ! test sample ULD interface.
96 :
97 : block
98 :
99 : type(lowDia_type), parameter :: subset = lowDia_type()
100 400 : dim = getChoice([1, 2])
101 400 : ndim = getUnifRand(1_IK, 5_IK)
102 6232 : cov_ref = getFilled(ZERO, ndim, ndim)
103 1978 : mean = getFilled(ZERO, ndim)
104 400 : if (dim == 2) then
105 3823 : sample = getUnifRand(ONE, TWO, ndim, nsam)
106 924 : meang = sample(:,1)
107 : else
108 3843 : sample = getUnifRand(ONE, TWO, nsam, ndim)
109 1054 : meang = sample(1,:)
110 : end if
111 2624 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
112 2624 : iweight = getUnifRand(1_IK, 9_IK, nsam)
113 2224 : rweisqs = sum(rweight**2)
114 :
115 : ! integer weighted
116 :
117 6232 : cov = getCov(sample, dim, iweight)
118 400 : call setCovMean(cov_ref, subset, mean, sample, dim, iweight, iweisum, meang)
119 400 : call setMatCopy(cov_ref, rdpack, cov_ref, rdpack, subset, transHerm)
120 400 : call report(__LINE__, SK_"integer-weighted")
121 :
122 11664 : cov = getCov(sample, dim, iweight, fweight_type())
123 5832 : cov_ref = cov_ref * getVarCorrection(real(iweisum, TKC))
124 400 : call report(__LINE__, SK_"integer-weighted")
125 :
126 : ! real weighted
127 :
128 6232 : cov = getCov(sample, dim, rweight)
129 400 : call setCovMean(cov_ref, subset, mean, sample, dim, rweight, rweisum, meang)
130 400 : call setMatCopy(cov_ref, rdpack, cov_ref, rdpack, subset, transHerm)
131 400 : call report(__LINE__, SK_"real-weighted")
132 :
133 11664 : cov = getCov(sample, dim, rweight, rweight_type())
134 5832 : cov_ref = cov_ref * getVarCorrection(rweisum, rweisqs)
135 400 : call report(__LINE__, SK_"real-weighted")
136 :
137 : ! unweighted
138 :
139 6232 : cov = getCov(sample, dim)
140 400 : call setCovMean(cov_ref, subset, mean, sample, dim, meang)
141 400 : call setMatCopy(cov_ref, rdpack, cov_ref, rdpack, subset, transHerm)
142 400 : call report(__LINE__, SK_"unweighted")
143 :
144 11664 : cov = getCov(sample, dim, rweight_type())
145 5832 : cov_ref = cov_ref * getVarCorrection(real(nsam, TKC))
146 400 : call report(__LINE__, SK_"unweighted")
147 :
148 : end block
149 :
150 : ! test sample XY interface.
151 :
152 8 : block
153 :
154 : type(lowDia_type), parameter :: subset = lowDia_type()
155 400 : ndim = 2_IK
156 1600 : mean = getFilled(ZERO, ndim)
157 3200 : cov_ref = getFilled(ZERO, ndim, ndim)
158 5248 : sample = getUnifRand(ONE, TWO, nsam, ndim)
159 2624 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
160 2624 : iweight = getUnifRand(1_IK, 9_IK, nsam)
161 2224 : rweisqs = sum(rweight**2)
162 1600 : meang = sample(1,:)
163 :
164 : ! integer weighted
165 :
166 3200 : cov = getCov(sample(:,1), sample(:,2), iweight)
167 400 : call setCovMean(cov_ref, mean, sample(:,1), sample(:,2), iweight, iweisum, meang)
168 400 : call report(__LINE__, SK_"integer-weighted")
169 :
170 5600 : cov = getCov(sample(:,1), sample(:,2), iweight, fweight_type())
171 2800 : cov_ref = cov_ref * getVarCorrection(real(iweisum, TKC))
172 400 : call report(__LINE__, SK_"integer-weighted")
173 :
174 : ! real weighted
175 :
176 3200 : cov = getCov(sample(:,1), sample(:,2), rweight)
177 400 : call setCovMean(cov_ref, mean, sample(:,1), sample(:,2), rweight, rweisum, meang)
178 400 : call report(__LINE__, SK_"real-weighted")
179 :
180 5600 : cov = getCov(sample(:,1), sample(:,2), rweight, rweight_type())
181 2800 : cov_ref = cov_ref * getVarCorrection(rweisum, rweisqs)
182 400 : call report(__LINE__, SK_"real-weighted")
183 :
184 : ! unweighted
185 :
186 3200 : cov = getCov(sample(:,1), sample(:,2))
187 400 : call setCovMean(cov_ref, mean, sample(:,1), sample(:,2), meang)
188 400 : call report(__LINE__, SK_"unweighted")
189 :
190 5600 : cov = getCov(sample(:,1), sample(:,2), rweight_type())
191 2800 : cov_ref = cov_ref * getVarCorrection(real(nsam, TKC))
192 400 : call report(__LINE__, SK_"unweighted")
193 :
194 : end block
195 :
196 : end do
197 :
198 : contains
199 :
200 5600 : subroutine setAssertedCov(cov)
201 : TYPE_OF_SAMPLE, intent(inout) :: cov(:,:)
202 98390 : cdiff = abs(cov - cov_ref)
203 64036 : assertion = assertion .and. all(cdiff < ctol)
204 5600 : end subroutine
205 :
206 800 : subroutine reportCFC(line)
207 : integer, intent(in) :: line
208 800 : call setAssertedCov(cov)
209 800 : if (test%traceable .and. .not. assertion) then
210 : ! LCOV_EXCL_START
211 : call test%disp%skip()
212 : call test%disp%show("ndim")
213 : call test%disp%show( ndim )
214 : call test%disp%show("std")
215 : call test%disp%show( std )
216 : call test%disp%show("cor")
217 : call test%disp%show( cor )
218 : call test%disp%show("cov")
219 : call test%disp%show( cov )
220 : call test%disp%show("cov_ref")
221 : call test%disp%show( cov_ref )
222 : call test%disp%show("cdiff")
223 : call test%disp%show( cdiff )
224 : call test%disp%skip()
225 : ! LCOV_EXCL_STOP
226 : end if
227 800 : call test%assert(assertion, SK_"The `cov` must be correctly computed from the specified `cor` and `std`.", int(line, IK))
228 800 : end subroutine
229 :
230 4800 : subroutine report(line, this)
231 : integer, intent(in) :: line
232 : character(*, SK), intent(in) :: this
233 4800 : call setAssertedCov(cov)
234 4800 : if (test%traceable .and. .not. assertion) then
235 : ! LCOV_EXCL_START
236 : call test%disp%skip()
237 : call test%disp%show("[ndim, nsam, dim, shape(sample, IK)]")
238 : call test%disp%show( [ndim, nsam, dim, shape(sample, IK)] )
239 : call test%disp%show("sample")
240 : call test%disp%show( sample )
241 : call test%disp%show("iweight")
242 : call test%disp%show( iweight )
243 : call test%disp%show("rweight")
244 : call test%disp%show( rweight )
245 : call test%disp%show("cov_ref")
246 : call test%disp%show( cov_ref )
247 : call test%disp%show("cov")
248 : call test%disp%show( cov )
249 : call test%disp%show("cdiff")
250 : call test%disp%show( cdiff )
251 : call test%disp%skip()
252 : ! LCOV_EXCL_STOP
253 : end if
254 4800 : call test%assert(assertion, SK_"The `cov` must be computed correctly for "//this//SK_" sample.", int(line, IK))
255 4800 : end subroutine
256 :
257 : !%%%%%%%%%%%%%
258 : #elif setCov_ENABLED
259 : !%%%%%%%%%%%%%
260 :
261 : real(TKC) :: rweisum
262 : integer(IK) :: iweisum
263 : real(TKC), allocatable :: rweight(:)
264 : integer(IK), allocatable :: iweight(:)
265 : integer(IK) :: itry, nsam, ndim, dim
266 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:)
267 : TYPE_OF_SAMPLE, allocatable :: cov(:,:), cor(:,:), cov_ref(:,:), cdiff(:,:)
268 : real(TKC), allocatable :: std(:)
269 8 : assertion = .true._LK
270 :
271 408 : do itry = 1, 50
272 :
273 : ! test CFC subsetr = lowDia.
274 :
275 : block
276 :
277 : type(lowDia_type), parameter :: subsetr = lowDia_type()
278 400 : ndim = getUnifRand(1_IK, 5_IK)
279 2020 : std = getUnifRand(1._TKC, 2._TKC, ndim)
280 6546 : cor = getUnifRand(-ONE, ONE, ndim, ndim)
281 400 : call setMatInit(cor, getSubSymm(subsetr), ZERO, RONE)
282 :
283 734 : cov_ref = getCFC(cor, subsetr, std)
284 :
285 6546 : cov = getFilled(ZERO, ndim, ndim)
286 400 : call setCov(cov, subsetr, cor, subsetr, std)
287 400 : call reportCFC(__LINE__, subsetr)
288 :
289 6546 : cov = getFilled(ZERO, ndim, ndim)
290 400 : call setCov(cov, getSubSymm(subsetr), cor, subsetr, std)
291 400 : call reportCFC(__LINE__, getSubSymm(subsetr))
292 :
293 : end block
294 :
295 : ! test CFC subsetr = uppDia.
296 :
297 : block
298 :
299 : type(uppDia_type), parameter :: subsetr = uppDia_type()
300 400 : ndim = getUnifRand(1_IK, 5_IK)
301 1998 : std = getUnifRand(1._TKC, 2._TKC, ndim)
302 6422 : cor = getUnifRand(-ONE, ONE, ndim, ndim)
303 400 : call setMatInit(cor, getSubSymm(subsetr), ZERO, RONE)
304 :
305 730 : cov_ref = getCFC(cor, subsetr, std)
306 :
307 6422 : cov = getFilled(ZERO, ndim, ndim)
308 400 : call setCov(cov, subsetr, cor, subsetr, std)
309 400 : call reportCFC(__LINE__, subsetr)
310 :
311 6422 : cov = getFilled(ZERO, ndim, ndim)
312 400 : call setCov(cov, getSubSymm(subsetr), cor, subsetr, std)
313 400 : call reportCFC(__LINE__, getSubSymm(subsetr))
314 :
315 : end block
316 :
317 400 : nsam = getUnifRand(2_IK, 7_IK)
318 :
319 : ! test sample lowDia interface.
320 :
321 : block
322 :
323 : type(lowDia_type), parameter :: subset = lowDia_type()
324 400 : dim = getChoice([1, 2])
325 400 : ndim = getUnifRand(1_IK, 5_IK)
326 400 : if (dim == 2) then
327 4375 : sample = getUnifRand(ONE, TWO, ndim, nsam)
328 : else
329 3645 : sample = getUnifRand(ONE, TWO, nsam, ndim)
330 : end if
331 2651 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
332 2651 : iweight = getUnifRand(1_IK, 9_IK, nsam)
333 2251 : iweisum = sum(iweight)
334 2251 : rweisum = sum(rweight)
335 :
336 : ! integer weighted
337 :
338 2005 : mean = getMean(sample, dim, iweight)
339 8355 : cov_ref = getCovRef(mean, sample, dim, real(iweight, TKC))
340 :
341 6504 : cov = getFilled(ZERO, ndim, ndim)
342 400 : call setCov(cov, subset, mean, sample, dim, iweight, iweisum)
343 400 : call report(__LINE__, subset, SK_"integer-weighted")
344 :
345 6504 : cov = getFilled(ZERO, ndim, ndim)
346 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim, iweight, iweisum)
347 400 : call report(__LINE__, subset, SK_"integer-weighted shifted")
348 :
349 : ! real weighted
350 :
351 2005 : mean = getMean(sample, dim, rweight)
352 6504 : cov_ref = getCovRef(mean, sample, dim, rweight)
353 :
354 6504 : cov = getFilled(ZERO, ndim, ndim)
355 400 : call setCov(cov, subset, mean, sample, dim, rweight, rweisum)
356 400 : call report(__LINE__, subset, SK_"real-weighted")
357 :
358 6504 : cov = getFilled(ZERO, ndim, ndim)
359 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim, rweight, rweisum)
360 400 : call report(__LINE__, subset, SK_"real-weighted shifted")
361 :
362 : ! unweighted
363 :
364 2005 : mean = getMean(sample, dim)
365 6504 : cov_ref = getCovRef(mean, sample, dim)
366 :
367 6504 : cov = getFilled(ZERO, ndim, ndim)
368 400 : call setCov(cov, subset, mean, sample, dim)
369 400 : call report(__LINE__, subset, SK_"unweighted")
370 :
371 6504 : cov = getFilled(ZERO, ndim, ndim)
372 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim)
373 400 : call report(__LINE__, subset, SK_"unweighted shifted")
374 :
375 : end block
376 :
377 : ! test sample uppDia interface.
378 :
379 : block
380 :
381 : type(uppDia_type), parameter :: subset = uppDia_type()
382 400 : dim = getChoice([1, 2])
383 400 : ndim = getUnifRand(1_IK, 5_IK)
384 400 : if (dim == 2) then
385 3647 : sample = getUnifRand(ONE, TWO, ndim, nsam)
386 : else
387 4196 : sample = getUnifRand(ONE, TWO, nsam, ndim)
388 : end if
389 2651 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
390 2651 : iweight = getUnifRand(1_IK, 9_IK, nsam)
391 2251 : iweisum = sum(iweight)
392 2251 : rweisum = sum(rweight)
393 :
394 : ! integer weighted
395 :
396 2005 : mean = getMean(sample, dim, iweight)
397 8349 : cov_ref = getCovRef(mean, sample, dim, real(iweight, TKC))
398 :
399 6498 : cov = getFilled(ZERO, ndim, ndim)
400 400 : call setCov(cov, subset, mean, sample, dim, iweight, iweisum)
401 400 : call report(__LINE__, subset, SK_"integer-weighted")
402 :
403 6498 : cov = getFilled(ZERO, ndim, ndim)
404 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim, iweight, iweisum)
405 400 : call report(__LINE__, subset, SK_"integer-weighted shifted")
406 :
407 : ! real weighted
408 :
409 2005 : mean = getMean(sample, dim, rweight)
410 6498 : cov_ref = getCovRef(mean, sample, dim, rweight)
411 :
412 6498 : cov = getFilled(ZERO, ndim, ndim)
413 400 : call setCov(cov, subset, mean, sample, dim, rweight, rweisum)
414 400 : call report(__LINE__, subset, SK_"real-weighted")
415 :
416 6498 : cov = getFilled(ZERO, ndim, ndim)
417 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim, rweight, rweisum)
418 400 : call report(__LINE__, subset, SK_"real-weighted shifted")
419 :
420 : ! unweighted
421 :
422 2005 : mean = getMean(sample, dim)
423 6498 : cov_ref = getCovRef(mean, sample, dim)
424 :
425 6498 : cov = getFilled(ZERO, ndim, ndim)
426 400 : call setCov(cov, subset, mean, sample, dim)
427 400 : call report(__LINE__, subset, SK_"unweighted")
428 :
429 6498 : cov = getFilled(ZERO, ndim, ndim)
430 1605 : call setCov(cov, subset, getShifted(sample, dim, -mean), dim)
431 400 : call report(__LINE__, subset, SK_"unweighted shifted")
432 :
433 : end block
434 :
435 : ! test sample XY interface.
436 :
437 8 : block
438 :
439 : type(lowDia_type), parameter :: subset = lowDia_type()
440 400 : dim = 1_IK
441 400 : ndim = 2_IK
442 1600 : mean = getFilled(ZERO, ndim)
443 3200 : cov_ref = getFilled(ZERO, ndim, ndim)
444 5302 : sample = getUnifRand(ONE, TWO, nsam, ndim)
445 2651 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
446 2651 : iweight = getUnifRand(1_IK, 9_IK, nsam)
447 2251 : iweisum = sum(iweight)
448 2251 : rweisum = sum(rweight)
449 :
450 : ! integer weighted
451 :
452 1600 : mean = getMean(sample, dim, iweight)
453 5051 : cov_ref = getCovRef(mean, sample, dim, real(iweight, TKC))
454 :
455 3200 : cov = getFilled(ZERO, ndim, ndim)
456 400 : call setCov(cov, mean, sample(:,1), sample(:,2), iweight, iweisum)
457 400 : call report(__LINE__, uppLowDia, SK_"integer-weighted")
458 :
459 3200 : cov = getFilled(ZERO, ndim, ndim)
460 4102 : call setCov(cov, sample(:,1) - mean(1), sample(:,2) - mean(2), iweight, iweisum)
461 400 : call report(__LINE__, uppLowDia, SK_"integer-weighted")
462 :
463 : ! real weighted
464 :
465 1600 : mean = getMean(sample, dim, rweight)
466 3200 : cov_ref = getCovRef(mean, sample, dim, rweight)
467 :
468 3200 : cov = getFilled(ZERO, ndim, ndim)
469 400 : call setCov(cov, mean, sample(:,1), sample(:,2), rweight, rweisum)
470 400 : call report(__LINE__, uppLowDia, SK_"real-weighted")
471 :
472 3200 : cov = getFilled(ZERO, ndim, ndim)
473 4102 : call setCov(cov, sample(:,1) - mean(1), sample(:,2) - mean(2), rweight, rweisum)
474 400 : call report(__LINE__, uppLowDia, SK_"real-weighted")
475 :
476 : ! unweighted
477 :
478 1600 : mean = getMean(sample, dim)
479 3200 : cov_ref = getCovRef(mean, sample, dim)
480 :
481 3200 : cov = getFilled(ZERO, ndim, ndim)
482 400 : call setCov(cov, mean, sample(:,1), sample(:,2))
483 400 : call report(__LINE__, uppLowDia, SK_"unweighted")
484 :
485 3200 : cov = getFilled(ZERO, ndim, ndim)
486 4102 : call setCov(cov, sample(:,1) - mean(1), sample(:,2) - mean(2))
487 400 : call report(__LINE__, uppLowDia, SK_"unweighted")
488 :
489 : end block
490 :
491 : end do
492 :
493 : contains
494 :
495 8800 : subroutine setAssertedCov(cov, subset)
496 : TYPE_OF_SAMPLE, intent(inout) :: cov(:,:)
497 : class(*), intent(in) :: subset
498 : select type(subset)
499 : type is (lowDia_type)
500 3200 : call setMatCopy(cov(1:,2:), rdpack, cov(2:,1:), rdpack, subset, transHerm)
501 : type is (uppDia_type)
502 3200 : call setMatCopy(cov(2:,1:), rdpack, cov(1:,2:), rdpack, subset, transHerm)
503 : type is (uppLowDia_type)
504 : continue
505 : class default
506 : error stop "Unrecognized subset." ! LCOV_EXCL_LINE
507 : end select
508 176412 : cdiff = abs(cov - cov_ref)
509 114348 : assertion = assertion .and. all(cdiff < ctol)
510 8800 : end subroutine
511 :
512 800 : pure function getCFC(cor, subsetr, std) result(cov)
513 : class(*), intent(in) :: subsetr
514 : real(TKC), intent(in) :: std(:)
515 : TYPE_OF_SAMPLE, intent(in) :: cor(:,:)
516 : TYPE_OF_SAMPLE :: cov(size(cor, 1, IK), size(cor, 2, IK))
517 : integer(IK) :: idim, jdim, ndim
518 800 : ndim = size(std, 1, IK)
519 3218 : do jdim = 1, ndim
520 8902 : do idim = jdim, ndim
521 8102 : if (same_type_as(subsetr, lowDia)) then
522 2873 : cov(idim, jdim) = cor(idim, jdim) * std(idim) * std(jdim)
523 2873 : cov(jdim, idim) = GET_CONJG(cov(idim, jdim))
524 2811 : elseif (same_type_as(subsetr, uppDia)) then
525 2811 : cov(jdim, idim) = cor(jdim, idim) * std(idim) * std(jdim)
526 2811 : cov(idim, jdim) = GET_CONJG(cov(jdim, idim))
527 : else
528 : error stop "Unrecognized subsetr." ! LCOV_EXCL_LINE
529 : end if
530 : end do
531 : end do
532 800 : end function
533 :
534 1600 : subroutine reportCFC(line, subset)
535 : class(*), intent(in) :: subset
536 : integer, intent(in) :: line
537 1600 : call setAssertedCov(cov, subset)
538 1600 : if (test%traceable .and. .not. assertion) then
539 : ! LCOV_EXCL_START
540 : call test%disp%skip()
541 : call test%disp%show("ndim")
542 : call test%disp%show( ndim )
543 : call test%disp%show("std")
544 : call test%disp%show( std )
545 : call test%disp%show("cor")
546 : call test%disp%show( cor )
547 : call test%disp%show("cov")
548 : call test%disp%show( cov )
549 : call test%disp%show("cov_ref")
550 : call test%disp%show( cov_ref )
551 : call test%disp%show("cdiff")
552 : call test%disp%show( cdiff )
553 : call test%disp%skip()
554 : ! LCOV_EXCL_STOP
555 : end if
556 1600 : call test%assert(assertion, SK_"The `cov` must be correctly computed from the specified `cor` and `std`.", int(line, IK))
557 1600 : end subroutine
558 :
559 3600 : PURE function getCovRef(mean, sample, dim, weight) result(cov)
560 : integer(IK), intent(in) :: dim
561 : real(TKC), intent(in), optional :: weight(:)
562 : TYPE_OF_SAMPLE, intent(in), contiguous :: sample(:,:), mean(:)
563 : TYPE_OF_SAMPLE :: cov(size(sample, 3 - dim, IK), size(sample, 3 - dim, IK))
564 7200 : TYPE_OF_SAMPLE :: sampleShifted(size(sample, 1, IK), size(sample, 2, IK))
565 : integer(IK) :: idim, jdim, ndim
566 : real(TKC) :: normfac
567 3600 : ndim = size(sample, 3 - dim, IK)
568 13230 : sampleShifted = getShifted(sample, dim, -mean)
569 3600 : if (present(weight)) then
570 13506 : normfac = 1._TKC / sum(weight)
571 2400 : if (dim == 2) then
572 3200 : do jdim = 1, ndim
573 11972 : do idim = 1, ndim
574 50912 : cov(idim, jdim) = dot_product(sampleShifted(jdim,:), sampleShifted(idim,:) * weight) * normfac
575 : end do
576 : end do
577 : else
578 5620 : do jdim = 1, ndim
579 18032 : do idim = 1, ndim
580 75670 : cov(idim, jdim) = dot_product(sampleShifted(:,jdim), sampleShifted(:,idim) * weight) * normfac
581 : end do
582 : end do
583 : end if
584 : else
585 1200 : normfac = 1._TKC / size(sample, dim, IK)
586 1200 : if (dim == 2) then
587 1600 : do jdim = 1, ndim
588 5986 : do idim = 1, ndim
589 25456 : cov(idim, jdim) = dot_product(sampleShifted(jdim,:), sampleShifted(idim,:)) * normfac
590 : end do
591 : end do
592 : else
593 2810 : do jdim = 1, ndim
594 9016 : do idim = 1, ndim
595 37835 : cov(idim, jdim) = dot_product(sampleShifted(:,jdim), sampleShifted(:,idim)) * normfac
596 : end do
597 : end do
598 : end if
599 : end if
600 3600 : end function
601 :
602 7200 : subroutine report(line, subset, this)
603 : integer, intent(in) :: line
604 : class(*), intent(in) :: subset
605 : character(*, SK), intent(in) :: this
606 7200 : call setAssertedCov(cov, subset)
607 7200 : if (test%traceable .and. .not. assertion) then
608 : ! LCOV_EXCL_START
609 : call test%disp%skip()
610 : call test%disp%show("[ndim, nsam, dim, shape(sample, IK)]")
611 : call test%disp%show( [ndim, nsam, dim, shape(sample, IK)] )
612 : call test%disp%show("sample")
613 : call test%disp%show( sample )
614 : call test%disp%show("iweight")
615 : call test%disp%show( iweight )
616 : call test%disp%show("rweight")
617 : call test%disp%show( rweight )
618 : call test%disp%show("cov_ref")
619 : call test%disp%show( cov_ref )
620 : call test%disp%show("cov")
621 : call test%disp%show( cov )
622 : call test%disp%show("cdiff")
623 : call test%disp%show( cdiff )
624 : call test%disp%skip()
625 : ! LCOV_EXCL_STOP
626 : end if
627 7200 : call test%assert(assertion, SK_"The `cov` must be computed correctly for "//this//SK_" sample.", int(line, IK))
628 7200 : end subroutine
629 :
630 : !%%%%%%%%%%%%%%%%%
631 : #elif setCovMean_ENABLED
632 : !%%%%%%%%%%%%%%%%%
633 :
634 : real(TKC) :: rweisum, rweisum_ref
635 : integer(IK) :: iweisum, iweisum_ref
636 : real(TKC), allocatable :: rweight(:)
637 : integer(IK), allocatable :: iweight(:)
638 : integer(IK) :: itry, nsam, ndim, dim
639 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:), meang(:), mean_ref(:), mdiff(:)
640 : TYPE_OF_SAMPLE, allocatable :: cov(:,:), cov_ref(:,:), cdiff(:,:)
641 8 : assertion = .true._LK
642 :
643 408 : do itry = 1, 50
644 :
645 400 : nsam = getUnifRand(2_IK, 7_IK)
646 :
647 : ! test sample lowDia interface.
648 :
649 : block
650 :
651 : type(lowDia_type), parameter :: subset = lowDia_type()
652 400 : ndim = getUnifRand(1_IK, 5_IK)
653 400 : dim = getChoice([1, 2])
654 400 : if (dim == 2) then
655 3988 : sample = getUnifRand(ONE, TWO, ndim, nsam)
656 1001 : meang = sample(:,1)
657 : else
658 3797 : sample = getUnifRand(ONE, TWO, nsam, ndim)
659 994 : meang = sample(1,:)
660 : end if
661 2625 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
662 2625 : iweight = getUnifRand(1_IK, 9_IK, nsam)
663 2225 : iweisum_ref = sum(iweight)
664 2225 : rweisum_ref = sum(rweight)
665 1995 : mean_ref = getFilled(ZERO, ndim)
666 6370 : cov_ref = getFilled(ZERO, ndim, ndim)
667 :
668 : ! integer weighted
669 :
670 1995 : mean_ref = getMean(sample, dim, iweight)
671 400 : call setCov(cov_ref, subset, mean_ref, sample, dim, iweight, iweisum_ref)
672 :
673 1995 : mean = getFilled(ZERO, ndim)
674 6370 : cov = getFilled(ZERO, ndim, ndim)
675 400 : call setCovMean(cov, subset, mean, sample, dim, iweight, iweisum, meang)
676 400 : call setAssertedSum(__LINE__, SK_"integer-weighted", real(iweisum, TKC), real(iweisum_ref, TKC))
677 400 : call setAssertedCov(__LINE__, SK_"integer-weighted", subset)
678 400 : call setAssertedAvg(__LINE__, SK_"integer-weighted")
679 :
680 : ! real weighted
681 :
682 1995 : mean_ref = getMean(sample, dim, rweight)
683 400 : call setCov(cov_ref, subset, mean_ref, sample, dim, rweight, rweisum_ref)
684 :
685 1995 : mean = getFilled(ZERO, ndim)
686 6370 : cov = getFilled(ZERO, ndim, ndim)
687 400 : call setCovMean(cov, subset, mean, sample, dim, rweight, rweisum, meang)
688 400 : call setAssertedSum(__LINE__, SK_"real-weighted", rweisum, rweisum_ref)
689 400 : call setAssertedCov(__LINE__, SK_"real-weighted", subset)
690 400 : call setAssertedAvg(__LINE__, SK_"real-weighted")
691 :
692 : ! unweighted
693 :
694 1995 : mean_ref = getMean(sample, dim)
695 400 : call setCov(cov_ref, subset, mean_ref, sample, dim)
696 :
697 1995 : mean = getFilled(ZERO, ndim)
698 6370 : cov = getFilled(ZERO, ndim, ndim)
699 400 : call setCovMean(cov, subset, mean, sample, dim, meang)
700 400 : call setAssertedCov(__LINE__, SK_"unweighted", subset)
701 400 : call setAssertedAvg(__LINE__, SK_"unweighted")
702 :
703 : end block
704 :
705 : ! test sample uppDia interface.
706 :
707 : block
708 :
709 : type(uppDia_type), parameter :: subset = uppDia_type()
710 400 : ndim = getUnifRand(1_IK, 5_IK)
711 400 : dim = getChoice([1, 2])
712 400 : if (dim == 2) then
713 4257 : sample = getUnifRand(ONE, TWO, ndim, nsam)
714 1043 : meang = sample(:,1)
715 : else
716 3337 : sample = getUnifRand(ONE, TWO, nsam, ndim)
717 903 : meang = sample(1,:)
718 : end if
719 2625 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
720 2625 : iweight = getUnifRand(1_IK, 9_IK, nsam)
721 2225 : iweisum_ref = sum(iweight)
722 2225 : rweisum_ref = sum(rweight)
723 1946 : mean_ref = getFilled(ZERO, ndim)
724 6046 : cov_ref = getFilled(ZERO, ndim, ndim)
725 :
726 : ! integer weighted
727 :
728 1946 : mean_ref = getMean(sample, dim, iweight)
729 400 : call setCov(cov_ref, subset, mean_ref, sample, dim, iweight, iweisum_ref)
730 :
731 1946 : mean = getFilled(ZERO, ndim)
732 6046 : cov = getFilled(ZERO, ndim, ndim)
733 400 : call setCovMean(cov, subset, mean, sample, dim, iweight, iweisum, meang)
734 400 : call setAssertedSum(__LINE__, SK_"integer-weighted", real(iweisum, TKC), real(iweisum_ref, TKC))
735 400 : call setAssertedCov(__LINE__, SK_"integer-weighted", subset)
736 400 : call setAssertedAvg(__LINE__, SK_"integer-weighted")
737 :
738 : ! real weighted
739 :
740 1946 : mean_ref = getMean(sample, dim, rweight)
741 400 : call setCov(cov_ref, subset, mean_ref, sample, dim, rweight, rweisum_ref)
742 :
743 1946 : mean = getFilled(ZERO, ndim)
744 6046 : cov = getFilled(ZERO, ndim, ndim)
745 400 : call setCovMean(cov, subset, mean, sample, dim, rweight, rweisum, meang)
746 400 : call setAssertedSum(__LINE__, SK_"real-weighted", rweisum, rweisum_ref)
747 400 : call setAssertedCov(__LINE__, SK_"real-weighted", subset)
748 400 : call setAssertedAvg(__LINE__, SK_"real-weighted")
749 :
750 : ! unweighted
751 :
752 1946 : mean_ref = getMean(sample, dim)
753 400 : call setCov(cov_ref, subset, mean_ref, sample, dim)
754 :
755 1946 : mean = getFilled(ZERO, ndim)
756 6046 : cov = getFilled(ZERO, ndim, ndim)
757 400 : call setCovMean(cov, subset, mean, sample, dim, meang)
758 400 : call setAssertedCov(__LINE__, SK_"unweighted", subset)
759 400 : call setAssertedAvg(__LINE__, SK_"unweighted")
760 :
761 : end block
762 :
763 : ! test sample XY interface.
764 :
765 8 : block
766 :
767 : type(uppLowDia_type), parameter :: subset = uppLowDia_type()
768 400 : dim = 1_IK
769 400 : ndim = 2_IK
770 5250 : sample = getUnifRand(ONE, TWO, nsam, ndim)
771 1600 : meang = sample(1,:)
772 2625 : rweight = getUnifRand(1._TKC, 9._TKC, nsam)
773 2625 : iweight = getUnifRand(1_IK, 9_IK, nsam)
774 2225 : iweisum_ref = sum(iweight)
775 2225 : rweisum_ref = sum(rweight)
776 1600 : mean_ref = getFilled(ZERO, ndim)
777 3200 : cov_ref = getFilled(ZERO, ndim, ndim)
778 :
779 : ! integer weighted
780 :
781 1600 : mean_ref = getMean(sample, dim, iweight)
782 400 : call setCov(cov_ref, mean_ref, sample(:,1), sample(:,2), iweight, iweisum_ref)
783 :
784 1600 : mean = getFilled(ZERO, ndim)
785 3200 : cov = getFilled(ZERO, ndim, ndim)
786 400 : call setCovMean(cov, mean, sample(:,1), sample(:,2), iweight, iweisum, meang)
787 400 : call setAssertedSum(__LINE__, SK_"integer-weighted", real(iweisum, TKC), real(iweisum_ref, TKC))
788 400 : call setAssertedCov(__LINE__, SK_"integer-weighted", subset)
789 400 : call setAssertedAvg(__LINE__, SK_"integer-weighted")
790 :
791 : ! real weighted
792 :
793 1600 : mean_ref = getMean(sample, dim, rweight)
794 400 : call setCov(cov_ref, mean_ref, sample(:,1), sample(:,2), rweight, rweisum_ref)
795 :
796 1600 : mean = getFilled(ZERO, ndim)
797 3200 : cov = getFilled(ZERO, ndim, ndim)
798 400 : call setCovMean(cov, mean, sample(:,1), sample(:,2), rweight, rweisum, meang)
799 400 : call setAssertedSum(__LINE__, SK_"real-weighted", rweisum, rweisum_ref)
800 400 : call setAssertedCov(__LINE__, SK_"real-weighted", subset)
801 400 : call setAssertedAvg(__LINE__, SK_"real-weighted")
802 :
803 : ! unweighted
804 :
805 1600 : mean_ref = getMean(sample, dim)
806 400 : call setCov(cov_ref, mean_ref, sample(:,1), sample(:,2))
807 :
808 1600 : mean = getFilled(ZERO, ndim)
809 3200 : cov = getFilled(ZERO, ndim, ndim)
810 400 : call setCovMean(cov, mean, sample(:,1), sample(:,2), meang)
811 400 : call setAssertedCov(__LINE__, SK_"unweighted", subset)
812 400 : call setAssertedAvg(__LINE__, SK_"unweighted")
813 :
814 : end block
815 :
816 : end do
817 :
818 : contains
819 :
820 2400 : subroutine setAssertedSum(line, this, weisum, weisum_ref)
821 : real(TKC), intent(in) :: weisum, weisum_ref
822 : character(*, SK), intent(in) :: this
823 : integer, intent(in) :: line
824 2400 : assertion = assertion .and. abs(weisum - weisum_ref) < rtol * weisum_ref
825 2400 : call report(line, SK_"The `weisum` must be computed correctly for "//this//SK_" sample.")
826 2400 : end subroutine
827 :
828 3600 : subroutine setAssertedAvg(line, this)
829 : character(*, SK), intent(in) :: this
830 : integer, intent(in) :: line
831 21381 : mdiff = abs(mean - mean_ref)
832 13023 : assertion = assertion .and. all(mdiff < ctol)
833 3600 : call report(line, SK_"The `mean` must be computed correctly for "//this//SK_" sample.")
834 3600 : end subroutine
835 :
836 3600 : subroutine setAssertedCov(line, this, subset)
837 : character(*, SK), intent(in) :: this
838 : class(*), intent(in) :: subset
839 : integer, intent(in) :: line
840 : select type(subset)
841 : type is (lowDia_type)
842 1200 : call setMatCopy(cov(1:,2:), rdpack, cov(2:,1:), rdpack, subset, transHerm)
843 1200 : call setMatCopy(cov_ref(1:,2:), rdpack, cov_ref(2:,1:), rdpack, subset, transHerm)
844 : type is (uppDia_type)
845 1200 : call setMatCopy(cov(2:,1:), rdpack, cov(1:,2:), rdpack, subset, transHerm)
846 1200 : call setMatCopy(cov_ref(2:,1:), rdpack, cov_ref(1:,2:), rdpack, subset, transHerm)
847 : type is (uppLowDia_type)
848 : continue
849 : class default
850 : error stop "Unrecognized subset." ! LCOV_EXCL_LINE
851 : end select
852 67182 : cdiff = abs(cov - cov_ref)
853 43248 : assertion = assertion .and. all(cdiff < ctol)
854 3600 : call report(line, SK_"The `cov` must be computed correctly for "//this//SK_" sample.")
855 3600 : end subroutine
856 :
857 9600 : subroutine report(line, msg)
858 : integer, intent(in) :: line
859 : character(*, SK), intent(in) :: msg
860 9600 : if (test%traceable .and. .not. assertion) then
861 : ! LCOV_EXCL_START
862 : call test%disp%skip()
863 : call test%disp%show("[ndim, nsam, dim, shape(sample, IK)]")
864 : call test%disp%show( [ndim, nsam, dim, shape(sample, IK)] )
865 : call test%disp%show("sample")
866 : call test%disp%show( sample )
867 : call test%disp%show("iweight")
868 : call test%disp%show( iweight )
869 : call test%disp%show("rweight")
870 : call test%disp%show( rweight )
871 : call test%disp%show("cov_ref")
872 : call test%disp%show( cov_ref )
873 : call test%disp%show("cov")
874 : call test%disp%show( cov )
875 : call test%disp%show("cdiff")
876 : call test%disp%show( cdiff )
877 : call test%disp%show("mean_ref")
878 : call test%disp%show( mean_ref )
879 : call test%disp%show("mean")
880 : call test%disp%show( mean )
881 : call test%disp%show("mdiff")
882 : call test%disp%show( mdiff )
883 : call test%disp%skip()
884 : ! LCOV_EXCL_STOP
885 : end if
886 9600 : call test%assert(assertion, msg, int(line, IK))
887 9600 : end subroutine
888 :
889 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
890 : #elif getCovMerged_ENABLED || setCovMerged_ENABLED
891 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
892 :
893 : real(TKC) :: fracA
894 : integer(IK) :: itry, ndim, dim
895 : integer(IK) :: nsam, nsamA, nsamB
896 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), meanDiff(:)
897 : TYPE_OF_SAMPLE, allocatable :: covA(:,:), covB(:,:), cov(:,:), cov_ref(:,:), cdiff(:,:)
898 16 : assertion = .true._LK
899 16 : dim = 2_IK
900 :
901 816 : do itry = 1, 50
902 :
903 800 : nsamA = getUnifRand(2_IK, 5_IK)
904 800 : nsamB = getUnifRand(2_IK, 5_IK)
905 800 : nsam = nsamA + nsamB
906 800 : fracA = real(nsamA, TKC) / real(nsam, TKC)
907 :
908 : ! test D2 interface.
909 :
910 16 : block
911 :
912 : type(uppDia_type), parameter :: subset = uppDia_type()
913 800 : ndim = getUnifRand(1_IK, 5_IK)
914 24095 : sample = getUnifRand(ONE, TWO, ndim, nsam)
915 :
916 12728 : cov_ref = getCov(sample, dim)
917 12728 : covA = getCov(sample(:,1:nsamA), dim)
918 12728 : covB = getCov(sample(:,nsamA+1:), dim)
919 7190 : meanDiff = getMean(sample(:,1:nsamA), dim) - getMean(sample(:,nsamA+1:), dim)
920 :
921 : #if getCovMerged_ENABLED
922 6302 : cov = getFilled(ZERO, ndim, ndim)
923 6302 : cov = getCovMerged(covB, covA, meanDiff, fracA)
924 8984 : cdiff = abs(cov - cov_ref)
925 5902 : assertion = assertion .and. all(cdiff < ctol)
926 400 : call report(__LINE__, SK_"The new `cov` must be computed correctly.")
927 : #elif setCovMerged_ENABLED
928 :
929 : ! new
930 :
931 6426 : cov = getFilled(ZERO, ndim, ndim)
932 400 : call setCovMerged(cov, covB, covA, meanDiff, fracA, subset)
933 400 : call setAssertedCov(__LINE__, subset, SK_"new")
934 :
935 6426 : cov = getFilled(ZERO, ndim, ndim)
936 400 : call setCovMerged(cov, covB, covA, meanDiff, fracA, getSubSymm(subset))
937 400 : call setAssertedCov(__LINE__, getSubSymm(subset), SK_"new")
938 :
939 : ! old
940 :
941 6426 : cov = covB
942 400 : call setCovMerged(cov, covA, meanDiff, fracA, subset)
943 400 : call setAssertedCov(__LINE__, subset, SK_"in-place")
944 :
945 6426 : cov = covB
946 400 : call setCovMerged(cov, covA, meanDiff, fracA, getSubSymm(subset))
947 400 : call setAssertedCov(__LINE__, getSubSymm(subset), SK_"in-place")
948 :
949 : #endif
950 :
951 : end block
952 :
953 : end do
954 :
955 : contains
956 :
957 : #if setCovMerged_ENABLED
958 1600 : subroutine setAssertedCov(line, subset, specific)
959 : character(*, SK), intent(in) :: specific
960 : class(*), intent(in) :: subset
961 : integer, intent(in) :: line
962 : select type (subset)
963 : type is (lowDia_type)
964 800 : call setMatCopy(cov, rdpack, cov, rdpack, subset, transHerm)
965 : type is (uppDia_type)
966 800 : call setMatCopy(cov, rdpack, cov, rdpack, subset, transHerm)
967 : class default
968 : error stop "Unrecognized subset." ! LCOV_EXCL_LINE
969 : end select
970 37088 : cdiff = abs(cov - cov_ref)
971 24104 : assertion = assertion .and. all(cdiff < ctol)
972 1600 : call report(line, SK_"The "//specific//SK_" `cov` must be computed correctly.")
973 1600 : end subroutine
974 : #endif
975 :
976 2000 : subroutine report(line, msg)
977 : integer, intent(in) :: line
978 : character(*, SK), intent(in) :: msg
979 2000 : if (test%traceable .and. .not. assertion) then
980 : ! LCOV_EXCL_START
981 : call test%disp%skip()
982 : call test%disp%show("[ndim, nsamA, nsamB, nsam, dim, shape(sample, IK)]")
983 : call test%disp%show( [ndim, nsamA, nsamB, nsam, dim, shape(sample, IK)] )
984 : call test%disp%show("sample")
985 : call test%disp%show( sample )
986 : call test%disp%show("meanDiff")
987 : call test%disp%show( meanDiff )
988 : call test%disp%show("covA")
989 : call test%disp%show( covA )
990 : call test%disp%show("covB")
991 : call test%disp%show( covB )
992 : call test%disp%show("cov_ref")
993 : call test%disp%show( cov_ref )
994 : call test%disp%show("cov")
995 : call test%disp%show( cov )
996 : call test%disp%show("cdiff")
997 : call test%disp%show( cdiff )
998 : call test%disp%skip()
999 : ! LCOV_EXCL_STOP
1000 : end if
1001 2000 : call test%assert(assertion, msg, int(line, IK))
1002 2000 : end subroutine
1003 :
1004 : !%%%%%%%%%%%%%%%%%%%%%%%
1005 : #elif setCovMeanMerged_ENABLED
1006 : !%%%%%%%%%%%%%%%%%%%%%%%
1007 :
1008 : real(TKC) :: fracA
1009 : integer(IK) :: itry, ndim, dim
1010 : integer(IK) :: nsam, nsamA, nsamB
1011 : TYPE_OF_SAMPLE, allocatable :: sample(:,:), mean(:), mean_ref(:), meanA(:), meanB(:), mdiff(:)
1012 : TYPE_OF_SAMPLE, allocatable :: covA_ref(:,:), covB_ref(:,:), covA(:,:), covB(:,:), cov(:,:), cov_ref(:,:), cdiff(:,:)
1013 8 : assertion = .true._LK
1014 8 : dim = 2_IK
1015 :
1016 408 : do itry = 1, 50
1017 :
1018 400 : nsamA = getUnifRand(2_IK, 5_IK)
1019 400 : nsamB = getUnifRand(2_IK, 5_IK)
1020 400 : nsam = nsamA + nsamB
1021 400 : fracA = real(nsamA, TKC) / real(nsam, TKC)
1022 :
1023 : ! test D2 interface.
1024 :
1025 8 : block
1026 :
1027 : type(uppDia_type), parameter :: subset = uppDia_type()
1028 400 : ndim = getUnifRand(1_IK, 5_IK)
1029 12203 : sample = getUnifRand(ONE, TWO, ndim, nsam)
1030 :
1031 6802 : cov_ref = getCov(sample, dim)
1032 2052 : mean_ref = getMean(sample, dim)
1033 6802 : covA_ref = getCov(sample(:,1:nsamA), dim)
1034 6802 : covB_ref = getCov(sample(:,nsamA+1:), dim)
1035 2052 : meanA = getMean(sample(:,1:nsamA), dim)
1036 2052 : meanB = getMean(sample(:,nsamA+1:), dim)
1037 :
1038 : ! new subset
1039 :
1040 6802 : covA = getFilled(ZERO, ndim, ndim); call setMatCopy(covA, rdpack, covA_ref, rdpack, subset)
1041 6802 : covB = getFilled(ZERO, ndim, ndim); call setMatCopy(covB, rdpack, covB_ref, rdpack, subset)
1042 6802 : cov = getFilled(ZERO, ndim, ndim)
1043 2052 : mean = getFilled(ZERO, ndim)
1044 :
1045 400 : call setCovMeanMerged(cov, mean, covB, meanB, covA, meanA, fracA, subset)
1046 400 : call setAssertedCov(__LINE__, subset, SK_"new")
1047 400 : call setAssertedAvg(__LINE__, SK_"new")
1048 :
1049 : ! new subset opposite
1050 :
1051 6802 : covA = getFilled(ZERO, ndim, ndim); call setMatCopy(covA, rdpack, covA_ref, rdpack, getSubSymm(subset))
1052 6802 : covB = getFilled(ZERO, ndim, ndim); call setMatCopy(covB, rdpack, covB_ref, rdpack, getSubSymm(subset))
1053 6802 : cov = getFilled(ZERO, ndim, ndim)
1054 2052 : mean = getFilled(ZERO, ndim)
1055 :
1056 400 : call setCovMeanMerged(cov, mean, covB, meanB, covA, meanA, fracA, getSubSymm(subset))
1057 400 : call setAssertedCov(__LINE__, getSubSymm(subset), SK_"new")
1058 400 : call setAssertedAvg(__LINE__, SK_"new")
1059 :
1060 : ! old subset
1061 :
1062 6802 : covA = getFilled(ZERO, ndim, ndim); call setMatCopy(covA, rdpack, covA_ref, rdpack, subset)
1063 6802 : covB = getFilled(ZERO, ndim, ndim); call setMatCopy(covB, rdpack, covB_ref, rdpack, subset)
1064 2052 : mean = meanB
1065 6802 : cov = covB
1066 :
1067 400 : call setCovMeanMerged(cov, mean, covA, meanA, fracA, subset)
1068 400 : call setAssertedCov(__LINE__, subset, SK_"in-place")
1069 400 : call setAssertedAvg(__LINE__, SK_"in-place")
1070 :
1071 : ! old subset opposite
1072 :
1073 6802 : covA = getFilled(ZERO, ndim, ndim); call setMatCopy(covA, rdpack, covA_ref, rdpack, getSubSymm(subset))
1074 6802 : covB = getFilled(ZERO, ndim, ndim); call setMatCopy(covB, rdpack, covB_ref, rdpack, getSubSymm(subset))
1075 2052 : mean = meanB
1076 6802 : cov = covB
1077 :
1078 400 : call setCovMeanMerged(cov, mean, covA, meanA, fracA, getSubSymm(subset))
1079 400 : call setAssertedCov(__LINE__, getSubSymm(subset), SK_"in-place")
1080 400 : call setAssertedAvg(__LINE__, SK_"in-place")
1081 :
1082 : end block
1083 :
1084 : end do
1085 :
1086 : contains
1087 :
1088 1600 : subroutine setAssertedAvg(line, specific)
1089 : character(*, SK), intent(in) :: specific
1090 : integer, intent(in) :: line
1091 10692 : mdiff = abs(mean - mean_ref)
1092 6608 : assertion = assertion .and. all(mdiff < ctol)
1093 1600 : call report(line, SK_"The "//specific//SK_" `meanMerged` must be computed correctly.")
1094 1600 : end subroutine
1095 :
1096 1600 : subroutine setAssertedCov(line, subset, specific)
1097 : character(*, SK), intent(in) :: specific
1098 : class(*), intent(in) :: subset
1099 : integer, intent(in) :: line
1100 : select type (subset)
1101 : type is (lowDia_type)
1102 800 : call setMatCopy(cov, rdpack, cov, rdpack, subset, transHerm)
1103 : type is (uppDia_type)
1104 800 : call setMatCopy(cov, rdpack, cov, rdpack, subset, transHerm)
1105 : class default
1106 : error stop "Unrecognized subset." ! LCOV_EXCL_LINE
1107 : end select
1108 39008 : cdiff = abs(cov - cov_ref)
1109 25608 : assertion = assertion .and. all(cdiff < ctol)
1110 1600 : call report(line, SK_"The "//specific//SK_" `cov` must be computed correctly.")
1111 1600 : end subroutine
1112 :
1113 3200 : subroutine report(line, msg)
1114 : integer, intent(in) :: line
1115 : character(*, SK), intent(in) :: msg
1116 3200 : if (test%traceable .and. .not. assertion) then
1117 : ! LCOV_EXCL_START
1118 : call test%disp%skip()
1119 : call test%disp%show("[ndim, nsamA, nsamB, nsam, dim, shape(sample, IK)]")
1120 : call test%disp%show( [ndim, nsamA, nsamB, nsam, dim, shape(sample, IK)] )
1121 : call test%disp%show("sample")
1122 : call test%disp%show( sample )
1123 : call test%disp%show("meanA")
1124 : call test%disp%show( meanA )
1125 : call test%disp%show("meanB")
1126 : call test%disp%show( meanB )
1127 : call test%disp%show("mean_ref")
1128 : call test%disp%show( mean_ref )
1129 : call test%disp%show("mean")
1130 : call test%disp%show( mean )
1131 : call test%disp%show("mdiff")
1132 : call test%disp%show( mdiff )
1133 : call test%disp%show("covA")
1134 : call test%disp%show( covA )
1135 : call test%disp%show("covB")
1136 : call test%disp%show( covB )
1137 : call test%disp%show("cov_ref")
1138 : call test%disp%show( cov_ref )
1139 : call test%disp%show("cov")
1140 : call test%disp%show( cov )
1141 : call test%disp%show("cdiff")
1142 : call test%disp%show( cdiff )
1143 : call test%disp%skip()
1144 : ! LCOV_EXCL_STOP
1145 : end if
1146 3200 : call test%assert(assertion, msg, int(line, IK))
1147 3200 : end subroutine
1148 :
1149 : #else
1150 : !%%%%%%%%%%%%%%%%%%%%%%%%
1151 : #error "Unrecognized interface."
1152 : !%%%%%%%%%%%%%%%%%%%%%%%%
1153 : #endif
1154 : #undef TYPE_OF_SAMPLE
1155 : #undef GET_CONJG
|
