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 the implementation details of the 2D routines under the generic interface [pm_sampleCor](@ref pm_sampleCor).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Saturday 4:40 PM, August 21, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the conjugation rule.
28 : #if CK_ENABLED
29 : #define TYPE_OF_SAMPLE complex(TKC)
30 : #define GET_CONJG(X)conjg(X)
31 : #define GET_RE(X)X%re
32 : #elif RK_ENABLED
33 : #define TYPE_OF_SAMPLE real(TKC)
34 : #define GET_CONJG(X)X
35 : #define GET_RE(X)X
36 : #elif getCor_ENABLED || setCor_ENABLED
37 : !elif !(SK_ENABLED || IK_ENABLED || LK_ENABLED || PSSK_ENABLED || BSSK_ENABLED)
38 : #error "Unrecognized interface."
39 : #endif
40 : ! Define weight type and kind and ZEROW.
41 : #if setRho_ENABLED && WTI_ENABLED
42 : #define TYPE_OF_WEIGHT integer(IK)
43 : #elif setRho_ENABLED && (WTR_ENABLED || WNO_ENABLED)
44 : #define TYPE_OF_WEIGHT real(RK)
45 : #endif
46 : ! Define the weight arguments.
47 : #if WNO_ENABLED
48 : #define WEIGHT_ARGS
49 : #elif WTI_ENABLED || WTR_ENABLED
50 : #define WEIGHT_ARGS, weight, weisum
51 : #elif (Prs_ENABLED || getRho_ENABLED || setRho_ENABLED)
52 : #error "Unrecognized interface."
53 : #endif
54 : !%%%%%%%%%%%%%%%%%%
55 : #if setCordance_ENABLED
56 : !%%%%%%%%%%%%%%%%%%
57 :
58 : ! Define the indexing rules.
59 : #if D0_ENABLED && SK_ENABLED
60 : #define GET_SIZE(ARRAY)len(ARRAY, IK)
61 : #define GETELL(ARRAY,I)ARRAY(I:I)
62 : #elif D1_ENABLED && (PSSK_ENABLED || BSSK_ENABLED)
63 : #define GET_SIZE(ARRAY)size(ARRAY, 1, IK)
64 : #define GETELL(ARRAY,I)ARRAY(I)%val
65 : #elif D1_ENABLED
66 : #define GET_SIZE(ARRAY)size(ARRAY, 1, IK)
67 : #define GETELL(ARRAY,I)ARRAY(I)
68 : #else
69 : #error "Unrecognized interface."
70 : #endif
71 : integer(IK) :: isam, jsam
72 372 : CHECK_ASSERTION(__LINE__, GET_SIZE(x) == GET_SIZE(y), SK_"@setCor(): The condition `size(x) == size(y)` must hold. size(x), size(y) = "//getStr([GET_SIZE(x), GET_SIZE(y)]))
73 124 : tiedx = 0_IK
74 124 : tiedy = 0_IK
75 : ! Set the output arguments.
76 : #if Sum_ENABLED
77 62 : cordance = 0_IK
78 : #define INCREMENT_CORDANCE cordance = cordance + 1_IK
79 : #define DECREMENT_CORDANCE cordance = cordance - 1_IK
80 : #elif All_ENABLED
81 62 : concordance = 0_IK
82 62 : discordance = 0_IK
83 : #define INCREMENT_CORDANCE concordance = concordance + 1_IK
84 : #define DECREMENT_CORDANCE discordance = discordance + 1_IK
85 : #else
86 : #error "Unrecognized interface."
87 : #endif
88 1252 : do isam = 2, GET_SIZE(x)
89 8856 : do jsam = 1, isam - 1
90 8732 : if (GETELL(x,isam) < GETELL(x,jsam)) then
91 3628 : if (GETELL(y,isam) < GETELL(y,jsam)) then
92 1772 : INCREMENT_CORDANCE
93 1856 : elseif (GETELL(y,jsam) < GETELL(y,isam)) then
94 1696 : DECREMENT_CORDANCE
95 : else
96 160 : tiedy = tiedy + 1_IK
97 : end if
98 3976 : elseif (GETELL(x,jsam) < GETELL(x,isam)) then
99 3666 : if (GETELL(y,isam) < GETELL(y,jsam)) then
100 1716 : INCREMENT_CORDANCE
101 1950 : elseif (GETELL(y,jsam) < GETELL(y,isam)) then
102 1808 : DECREMENT_CORDANCE
103 : else
104 142 : tiedy = tiedy + 1_IK
105 : end if
106 : else
107 310 : tiedx = tiedx + 1_IK
108 310 : if (GETELL(y,isam) == GETELL(y,jsam)) tiedy = tiedy + 1_IK
109 : end if
110 : end do
111 : end do
112 : #undef INCREMENT_CORDANCE
113 : #undef DECREMENT_CORDANCE
114 : #undef GET_SIZE
115 : #undef GETELL
116 :
117 : ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118 : !#elif getCorMerged_ENABLED && New_ENABLED
119 : ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120 : !
121 : ! integer(IK) :: idim
122 : ! call setCorMerged(corMerged, corB, corA, meanDiff, fracA, uppDia)
123 : ! ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
124 : ! !do concurrent(idim = 2 : size(corMerged, 1, IK))
125 : ! do idim = 2, size(corMerged, 1, IK)
126 : ! corMerged(idim, 1 : idim - 1) = corMerged(1 : idim - 1, idim)
127 : ! end do
128 : !
129 : ! !%%%%%%%%%%%%%%%%%%%
130 : !#elif setCorMerged_ENABLED
131 : ! !%%%%%%%%%%%%%%%%%%%
132 : !
133 : ! integer(IK) :: idim, jdim, ndim
134 : ! real(TKC) :: fracB, fracAB, stdinv(size(meanDiff, 1, IK))
135 : ! ! Define the loop ranges for different subsets.
136 : !#if UXD_ENABLED || UXX_ENABLED
137 : !#define GET_RANGE(I,J,K)I, J - I
138 : !#elif XLD_ENABLED || XLX_ENABLED
139 : !#define GET_RANGE(I,J,K)J + I, K
140 : !#else
141 : !#error "Unrecognized interface."
142 : !#endif
143 : ! ! Define the output value for setCorMerged.
144 : !#if Old_ENABLED
145 : !#define corMerged corB
146 : !#elif New_ENABLED
147 : ! CHECK_ASSERTION(__LINE__, all(shape(corMerged, IK) == shape(corA, IK)), SK_"@setCorMerged(): The condition `all(shape(corMerged) == shape(corA))` must hold. shape(corMerged), shape(corA) = "//getStr([shape(corMerged, IK), shape(corA, IK)]))
148 : !#else
149 : !#error "Unrecognized interface."
150 : !#endif
151 : ! CHECK_ASSERTION(__LINE__, 0._TKC < fracA .and. fracA < 1._TKC, SK_"@setCorMerged(): The condition `0 < fracA .and. fracA < 1` must hold. fracA = "//getStr(fracA))
152 : ! CHECK_ASSERTION(__LINE__, all(shape(corB, IK) == shape(corA, IK)), SK_"@setCorMerged(): The condition `all(shape(corB, IK) == shape(corA, IK))` must hold. shape(corB), shape(corA) = "//getStr([shape(corB, IK), shape(corA, IK)]))
153 : ! CHECK_ASSERTION(__LINE__, all(size(meanDiff, 1, IK) == shape(corA, IK)), SK_"@setCorMerged(): The condition `all(size(meanDiff) == shape(corA))` must hold. size(meanDiff), shape(corA) = "//getStr([size(meanDiff, 1, IK), shape(corA, IK)]))
154 : ! fracB = 1._TKC - fracA
155 : ! fracAB = fracA * fracB
156 : ! ndim = size(corMerged, 1, IK)
157 : ! do jdim = 1, size(corA, 1, IK)
158 : ! stdinv(jdim) = 1._TKC / sqrt(1._TKC + fracAB * meanDiff(jdim)**2)
159 : !#if XLD_ENABLED
160 : ! corMerged(jdim, jdim) = 1._TKC
161 : !#endif
162 : ! do idim = GET_RANGE(1,jdim,ndim)
163 : ! corMerged(idim, jdim) = fracB * corB(idim, jdim) + fracA * corA(idim, jdim) + fracAB * meanDiff(idim) * meanDiff(jdim)
164 : ! end do
165 : !#if UXD_ENABLED
166 : ! corMerged(jdim, jdim) = 1._TKC
167 : !#endif
168 : ! end do
169 : ! do jdim = 1, ndim
170 : ! do idim = GET_RANGE(1,jdim,ndim)
171 : ! corMerged(idim, jdim) = corMerged(idim, jdim) * stdinv(idim) * stdinv(jdim)
172 : ! end do
173 : ! end do
174 : !#undef GET_RANGE
175 : !#undef corMerged
176 : !
177 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
178 : #elif getCor_ENABLED && CFC_ENABLED && RULD_ENABLED
179 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 :
181 : integer(IK) :: idim
182 : #if VUXX_ENABLED || VXLX_ENABLED
183 820 : call setCor(cor, uppDia, cov, subsetv, stdinv)
184 : #elif VUXD_ENABLED || VXLD_ENABLED
185 875 : call setCor(cor, uppDia, cov, subsetv)
186 : #else
187 : #error "Unrecognized interface."
188 : #endif
189 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
190 : !do concurrent(idim = 2 : size(cor, 1, IK))
191 4993 : do idim = 2, size(cor, 1, IK)
192 11741 : cor(idim, 1 : idim - 1) = GET_CONJG(cor(1 : idim - 1, idim))
193 : end do
194 :
195 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196 : #elif getCor_ENABLED && (XY_ENABLED || ULD_ENABLED)
197 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
198 :
199 : #if WTI_ENABLED
200 : integer(IK) :: weisum
201 : #elif WTR_ENABLED
202 : real(TKC) :: weisum
203 : #elif !WNO_ENABLED
204 : #error "Unrecognized interface."
205 : #endif
206 : #if XY_ENABLED
207 : TYPE_OF_SAMPLE :: mean(2)
208 4803 : call setMean(mean, x, y WEIGHT_ARGS)
209 4803 : call setCor(cor, mean, x, y WEIGHT_ARGS)
210 : #elif ULD_ENABLED
211 12312 : TYPE_OF_SAMPLE :: mean(size(sample, 3 - dim, IK))
212 6156 : call setMean(mean, sample, dim WEIGHT_ARGS)
213 6156 : call setCor(cor, uppDia, mean, sample, dim WEIGHT_ARGS)
214 : #else
215 : #error "Unrecognized interface."
216 : #endif
217 : #if ULD_ENABLED
218 : block
219 : integer(IK) :: ndim, idim
220 6156 : ndim = size(cor, 1, IK)
221 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
222 : !do concurrent(idim = 1 : ndim)
223 24669 : do idim = 1, ndim
224 49635 : cor(idim, 1 : idim - 1) = GET_CONJG(cor(1 : idim - 1, idim))
225 : end do
226 : end block
227 : #endif
228 :
229 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
230 : #elif setCor_ENABLED && CFC_ENABLED
231 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232 :
233 : integer(IK) :: idim, jdim, ndim
234 : ! Define the indexing rules for `cor`.
235 : #if RUXX_ENABLED || RUXD_ENABLED
236 : #define RINDEX(I,J) I,J
237 : #elif RXLX_ENABLED || RXLD_ENABLED
238 : #define RINDEX(I,J) J,I
239 : #else
240 : #error "Unrecognized interface."
241 : #endif
242 : ! Define the indexing rules for `cov`.
243 : #if VUXX_ENABLED || VUXD_ENABLED
244 : #define VINDEX(I,J) I,J
245 : #elif VXLX_ENABLED || VXLD_ENABLED
246 : #define VINDEX(I,J) J,I
247 : #else
248 : #error "Unrecognized interface."
249 : #endif
250 : ! Define the conjugation rule.
251 : #if (VUXD_ENABLED && (RUXD_ENABLED || RUXX_ENABLED)) || (VUXX_ENABLED && (RUXD_ENABLED || RUXX_ENABLED)) || \
252 : (VXLD_ENABLED && (RXLD_ENABLED || RXLX_ENABLED)) || (VXLX_ENABLED && (RXLD_ENABLED || RXLX_ENABLED))
253 : #define CONJUGATE(X)X
254 : #elif (VUXD_ENABLED && (RXLD_ENABLED || RXLX_ENABLED)) || (VUXX_ENABLED && (RXLD_ENABLED || RXLX_ENABLED)) || \
255 : (VXLD_ENABLED && (RUXD_ENABLED || RUXX_ENABLED)) || (VXLX_ENABLED && (RUXD_ENABLED || RUXX_ENABLED))
256 : #define CONJUGATE(X)GET_CONJG(X)
257 : #else
258 : #error "Unrecognized interface."
259 : #endif
260 : ! Define the diagonal elements of `cor` if needed.
261 : #if (RUXD_ENABLED || RXLD_ENABLED) && (VUXD_ENABLED || VXLD_ENABLED)
262 : #define SET_CORDIA_BEG(DIM) cor(DIM, DIM) = 1._TKC
263 : #define SET_CORDIA_END(DIM)
264 : #elif (RUXD_ENABLED || RXLD_ENABLED) && (VUXX_ENABLED || VXLX_ENABLED)
265 : #define SET_CORDIA_BEG(DIM)
266 : #define SET_CORDIA_END(DIM) cor(DIM, DIM) = 1._TKC
267 : #elif RUXX_ENABLED || RXLX_ENABLED
268 : #define SET_CORDIA_BEG(DIM)
269 : #define SET_CORDIA_END(DIM)
270 : #else
271 : #error "Unrecognized interface."
272 : #endif
273 : ! Define the inverse standard deviations if needed.
274 : #if VUXD_ENABLED || VXLD_ENABLED
275 : #define SET_STDINV(DIM) stdinv(DIM) = 1._TKC / sqrt(GET_RE(cov(DIM, DIM)));
276 : #define SET_RANGE(START,STOP) STOP,START, -1
277 8310 : real(TKC) :: stdinv(size(cor, 1, IK))
278 : !do concurrent(idim = 1 : size(stdinv, 1, IK); stdinv(idim) = 1._TKC / sqrt(cov(idim, idim)); end do
279 : #elif VUXX_ENABLED || VXLX_ENABLED
280 : #define SET_RANGE(START,STOP) START,STOP
281 : #define SET_STDINV(DIM)
282 12300 : CHECK_ASSERTION(__LINE__, size(cov, 1, IK) == size(stdinv, 1, IK), SK_"@setCor(): The condition `size(cov, 1) == size(stdinv)` must hold. size(cov, 1), size(stdinv) = "//getStr([size(cov, 1, IK), size(stdinv, 1, IK)]))
283 16213 : CHECK_ASSERTION(__LINE__, all(0._TKC < stdinv), SK_"@setCor(): The condition `all(0. < stdinv)` must hold. stdinv = "//getStr(stdinv))
284 : #else
285 : #error "Unrecognized interface."
286 : #endif
287 4155 : ndim = size(cov, 1, IK)
288 8255 : CHECK_ASSERTION(__LINE__, ndim == size(cov, 2, IK), SK_"@setCor(): The condition `size(cov, 1) == size(cov, 2)` must hold. shape(cov) = "//getStr(shape(cov, IK)))
289 66040 : CHECK_ASSERTION(__LINE__, all(ndim == shape(cor, IK)), SK_"@setCor(): The condition `all(size(cov, 1) == shape(cor))` must hold. size(cov, 1), shape(cor) = "//getStr([ndim, shape(cor, IK)]))
290 32648 : do jdim = 1, ndim
291 12280 : SET_STDINV(jdim)
292 7364 : SET_CORDIA_BEG(jdim)
293 62658 : do idim = SET_RANGE(1, jdim - 1)
294 56863 : cor(RINDEX(idim, jdim)) = stdinv(idim) * stdinv(jdim) * CONJUGATE(cov(VINDEX(idim, jdim)))
295 : end do
296 9657 : SET_CORDIA_END(jdim)
297 : end do
298 :
299 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
300 : #elif setCor_ENABLED && Prs_ENABLED && XY_ENABLED
301 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
302 :
303 : TYPE_OF_SAMPLE :: cov(2,2)
304 : #if Org_ENABLED
305 1203 : call setCov(cov, x, y WEIGHT_ARGS)
306 : #elif Avg_ENABLED
307 10815 : call setCov(cov, mean, x, y WEIGHT_ARGS)
308 : #else
309 : #error "Unrecognized interface."
310 : #endif
311 12018 : cor = cov(1,2) / sqrt(GET_RE(cov(1,1)) * GET_RE(cov(2,2)))
312 :
313 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
314 : #elif setCor_ENABLED && Prs_ENABLED && (UXD_ENABLED || XLD_ENABLED)
315 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
316 :
317 : integer(IK) :: idim, jdim, ndim
318 34306 : real(TKC) :: stdinv(size(cor, 1, IK))
319 : #if Org_ENABLED
320 2413 : call setCov(cor, subset, sample, dim WEIGHT_ARGS)
321 : #elif Avg_ENABLED
322 14740 : call setCov(cor, subset, mean, sample, dim WEIGHT_ARGS)
323 : #else
324 : #error "Unrecognized interface."
325 : #endif
326 : ndim = size(cor, 1, IK)
327 68763 : do idim = 1, ndim
328 51610 : stdinv(idim) = 1._TKC / sqrt(GET_RE(cor(idim, idim)))
329 68763 : cor(idim, idim) = 1._TKC
330 : end do
331 : #if UXD_ENABLED
332 : #define THIS_RANGE 1, jdim - 1
333 : #elif XLD_ENABLED
334 : #define THIS_RANGE jdim + 1, ndim
335 : #else
336 : #error "Unrecognized interface."
337 : #endif
338 68763 : do jdim = 1, ndim
339 138098 : do idim = THIS_RANGE
340 120945 : cor(idim, jdim) = cor(idim, jdim) * stdinv(idim) * stdinv(jdim)
341 : end do
342 : end do
343 : #undef THIS_RANGE
344 :
345 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
346 : #elif getRho_ENABLED && (XY_D0_ENABLED || XY_D1_ENABLED)
347 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
348 :
349 : #if XY_D0_ENABLED
350 302 : real(RK) :: frankx(len(x, IK)), franky(len(x, IK))
351 : #elif XY_D1_ENABLED
352 3310 : real(RK) :: frankx(size(x, 1, IK)), franky(size(x, 1, IK))
353 : #else
354 : #error "Unrecognized interface."
355 : #endif
356 : #if WNO_ENABLED
357 604 : call setRho(rho, frankx, franky, x, y)
358 : #elif WTI_ENABLED || WTR_ENABLED
359 1202 : call setRho(rho, frankx, franky, x, y, weight)
360 : #else
361 : #error "Unrecognized interface."
362 : #endif
363 :
364 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
365 : #elif getRho_ENABLED && ULD_ENABLED
366 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%
367 :
368 : integer(IK) :: idim, ndim
369 1660 : real(RK) :: frank(size(sample, 1, IK), size(sample, 2, IK))
370 : #if WNO_ENABLED
371 556 : call setRho(rho, uppDia, frank, sample, dim)
372 : #elif WTI_ENABLED || WTR_ENABLED
373 1104 : call setRho(rho, uppDia, frank, sample, dim, weight)
374 : #else
375 : #error "Unrecognized interface."
376 : #endif
377 1660 : ndim = size(rho, 1, IK)
378 : ! The `do concurrent` version is problematic for OpenMP builds of the ParaMonte library with the Intel ifort 2021.8 on heap, yielding segmentation fault.
379 : !do concurrent(idim = 1 : ndim)
380 6819 : do idim = 1, ndim
381 13885 : rho(idim, 1 : idim - 1) = rho(1 : idim - 1, idim)
382 : end do
383 :
384 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
385 : #elif setRho_ENABLED && (XY_D0_ENABLED || XY_D1_ENABLED)
386 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
387 :
388 : real(RK) :: mean(2)
389 : #if WTI_ENABLED || WTR_ENABLED
390 : TYPE_OF_WEIGHT :: weisum
391 : #elif !WNO_ENABLED
392 : #error "Unrecognized interface."
393 : #endif
394 : #if XY_D0_ENABLED
395 903 : CHECK_ASSERTION(__LINE__, size(frankx, 1, IK) == len(x, IK), SK_"@setRho(): The condition `size(frankx) == len(x)` must hold. size(frankx), size(x) = "//getStr([size(frankx, 1, IK), len(x, IK)]))
396 903 : CHECK_ASSERTION(__LINE__, size(franky, 1, IK) == len(y, IK), SK_"@setRho(): The condition `size(franky) == len(y)` must hold. size(franky), size(y) = "//getStr([size(franky, 1, IK), len(y, IK)]))
397 : #elif XY_D1_ENABLED
398 9924 : CHECK_ASSERTION(__LINE__, size(frankx, 1, IK) == size(x, 1, IK), SK_"@setRho(): The condition `size(frankx) == size(x)` must hold. size(frankx), size(x) = "//getStr([size(frankx, 1, IK), size(x, 1, IK)]))
399 9924 : CHECK_ASSERTION(__LINE__, size(franky, 1, IK) == size(y, 1, IK), SK_"@setRho(): The condition `size(franky) == size(y)` must hold. size(franky), size(y) = "//getStr([size(franky, 1, IK), size(y, 1, IK)]))
400 : #else
401 : #error "Unrecognized interface."
402 : #endif
403 3609 : call setRankFractional(frankx, x)
404 3609 : call setRankFractional(franky, y)
405 3609 : call setMean(mean, frankx, franky WEIGHT_ARGS)
406 3609 : call setCor(rho, mean, frankx, franky WEIGHT_ARGS)
407 :
408 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
409 : #elif setRho_ENABLED && (UXD_ENABLED || XLD_ENABLED)
410 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
411 :
412 : #if WTI_ENABLED || WTR_ENABLED
413 : TYPE_OF_WEIGHT :: weisum
414 : #elif !WNO_ENABLED
415 : #error "Unrecognized interface."
416 : #endif
417 : integer(IK) :: idim, ndim, nsam
418 9944 : real(RK) :: mean(size(sample, 3 - dim, IK))
419 54692 : CHECK_ASSERTION(__LINE__, all(shape(frank, IK) == shape(sample, IK)), SK_"@setRho(): The condition `all(shape(frank) == shape(sample))` must hold. shape(frank), shape(sample) = "//getStr([shape(frank, IK), shape(sample, IK)]))
420 4972 : ndim = size(rho, 3 - dim, IK)
421 4972 : nsam = size(sample, dim, IK)
422 4972 : if (dim == 1_IK) then
423 9981 : do idim = 1, ndim
424 9981 : call setRankFractional(frank(1:nsam, idim), sample(1:nsam, idim))
425 : end do
426 : else
427 10020 : do idim = 1, ndim
428 385330 : call setRankFractional(frank(idim, 1:nsam), sample(idim, 1:nsam))
429 : end do
430 : end if
431 4972 : call setMean(mean, frank, dim WEIGHT_ARGS)
432 4972 : call setCor(rho, subset, mean, frank, dim WEIGHT_ARGS)
433 :
434 : #else
435 : !%%%%%%%%%%%%%%%%%%%%%%%%
436 : #error "Unrecognized interface."
437 : !%%%%%%%%%%%%%%%%%%%%%%%%
438 : #endif
439 : #undef TYPE_OF_SAMPLE
440 : #undef TYPE_OF_WEIGHT
441 : #undef SET_CORDIA_BEG
442 : #undef SET_CORDIA_END
443 : #undef WEIGHT_ARGS
444 : #undef SET_STDINV
445 : #undef CONJUGATE
446 : #undef GET_CONJG
447 : #undef SET_RANGE
448 : #undef GET_RE
449 : #undef RINDEX
450 : #undef VINDEX
|
