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 routines of [pm_sampleVar](@ref pm_sampleVar).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Saturday 4:40 PM, August 21, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Define the looping ranges for the corresponding matrix subsets.
28 : #if (setCov_ENABLED || setCovMean_ENABLED || setCovMerged_ENABLED || setCovMeanMerged_ENABLED) && (XLD_ENABLED || XLD_XLD_ENABLED || XLD_UXD_ENABLED)
29 : ! Start from the beginning of the lower-triangle.
30 : #define OFF_RANGE(I,J,OFFSET)I + OFFSET, J
31 : #define ROW_RANGE(I,J,K)J + 1_IK, K
32 : #define COL_RANGE(I,J)I, J
33 : #define FIRST 1
34 : #elif (setCov_ENABLED || setCovMean_ENABLED || setCovMerged_ENABLED || setCovMeanMerged_ENABLED) && (UXD_ENABLED || UXD_UXD_ENABLED || UXD_XLD_ENABLED)
35 : ! Start from the end of the upper-triangle.
36 : #define OFF_RANGE(I,J,OFFSET)J - OFFSET, I, -1_IK
37 : #define ROW_RANGE(I,J,K)J - 1_IK, I, -1_IK
38 : #define COL_RANGE(I,J)J, I, -1_IK
39 : #define FIRST ndim
40 : #elif setCovMerged_ENABLED || (setCov_ENABLED && !(XY_ENABLED || CorStd_ENABLED))
41 : #error "Unrecognized interface."
42 : #endif
43 : ! Set the conjugation rule.
44 : #if CK_ENABLED
45 : #define GET_RE(X)X%re
46 : #define SET_CONJG(X)X = conjg(X)
47 : #define GET_CONJG(X)conjg(X)
48 : #define TYPE_OF_SAMPLE complex(TKC)
49 : #define GET_ABSQ(X)(real(X)**2 + aimag(X)**2)
50 : #define GET_PROD(X,Y)(X%re * Y%re + X%im * Y%im)
51 : #elif RK_ENABLED
52 : #define GET_RE(X)X
53 : #define SET_CONJG(X)
54 : #define GET_CONJG(X)X
55 : #define TYPE_OF_SAMPLE real(TKC)
56 : #define GET_PROD(X,Y)X * Y
57 : #define GET_ABSQ(X)X**2
58 : #else
59 : #error "Unrecognized interface."
60 : #endif
61 : ! Set the shifting rule.
62 : #if setCov_ENABLED && Org_ENABLED
63 : #define GET_SHIFTED(X,Y)X
64 : #elif setCov_ENABLED && Avg_ENABLED
65 : #define GET_SHIFTED(X,Y)(X - Y)
66 : #elif setCov_ENABLED && !CorStd_ENABLED
67 : #error "Unrecognized interface."
68 : #endif
69 : ! Set the weighting rule.
70 : #if (setCov_ENABLED || setCovMean_ENABLED) && WNO_ENABLED
71 : #define GET_WEIGHTED(X,Y)X
72 : #elif (setCov_ENABLED || setCovMean_ENABLED) && (WTI_ENABLED || WTR_ENABLED)
73 : #define GET_WEIGHTED(X,Y)X * Y
74 : #elif setCovMean_ENABLED || (setCov_ENABLED && !CorStd_ENABLED)
75 : #error "Unrecognized interface."
76 : #endif
77 : ! Define weight type and kind and ZEROW.
78 : #if WTI_ENABLED
79 : #define TYPE_OF_WEIGHT integer(IK)
80 : #define ZEROW 0_IK
81 : #elif WTR_ENABLED || WNO_ENABLED
82 : #define TYPE_OF_WEIGHT real(TKC)
83 : #define ZEROW 0._TKC
84 : #elif (getCov_ENABLED || setCov_ENABLED) && !(WNO_ENABLED || CorStd_ENABLED)
85 : #error "Unrecognized interface."
86 : #endif
87 : ! Define runtime sanity checks.
88 : #define CHECK_VAL_NSAM(PROC,DIM) \
89 : CHECK_ASSERTION(__LINE__, 1 < size(sample, DIM, IK), \
90 : PROC//SK_": The condition `1 < size(sample, dim)` must hold. dim, shape(sample) = "//getStr([DIM, shape(sample, IK)]))
91 : #define CHECK_VAL_NDIM(PROC,DIM) \
92 : CHECK_ASSERTION(__LINE__, 0 < size(sample, DIM, IK), \
93 : PROC//SK_": The condition `0 < size(sample, dim)` must hold. dim, shape(sample) = "//getStr([DIM, shape(sample, IK)]))
94 : #define CHECK_VAL_DIM(PROC) \
95 : CHECK_ASSERTION(__LINE__, 1 <= dim .and. dim <= rank(sample), \
96 : PROC//SK_": The condition `1 <= dim .and. dim <= rank(sample)` must hold. dim, rank(sample) = "\
97 : //getStr([integer(IK) :: dim, rank(sample)]))
98 : #define CHECK_SUM_WEI(PROC) \
99 : CHECK_ASSERTION(__LINE__, ZEROW < sum(weight), \
100 : PROC//SK_": The condition `0 < sum(weight)` must hold. weight = "//getStr(weight))
101 : #define CHECK_VAL_WEI(PROC) \
102 : CHECK_ASSERTION(__LINE__, all(0._TKC <= weight), \
103 : PROC//SK_": The condition `all(0. <= weight)` must hold. weight = "//getStr(weight))
104 : #define CHECK_SHAPE_COV(PROC) \
105 : CHECK_ASSERTION(__LINE__, all(size(sample, 3 - dim, IK) == shape(cov, IK)), \
106 : PROC//SK_": The condition `all(size(sample, 3 - dim) == shape(cov))` must hold. dim, size(sample, 3 - dim), shape(cov) = "\
107 : //getStr([dim, size(sample, 3 - dim, IK), shape(cov, IK)]))
108 : #define CHECK_VAL_MEANG(PROC,dim) \
109 : CHECK_ASSERTION(__LINE__, all([minval(sample, dim) <= meang .and. meang <= maxval(sample, dim)]), \
110 : PROC//SK_": The condition `all([minval(sample, dim) <= meang .and. meang <= maxval(sample, dim)])` must hold. dim, minval(sample, dim), meang, maxval(sample, dim) = "//\
111 : getStr(dim)//SK_"; "//getStr(minval(sample, dim))//SK_"; "//getStr(meang)//SK_"; "//getStr(maxval(sample, dim)))
112 : #define CHECK_LEN_MEANG(PROC) \
113 : CHECK_ASSERTION(__LINE__, size(sample, 3 - dim, IK) == size(meang, 1, IK), \
114 : PROC//SK_": The condition `size(sample, 3 - dim) == size(meang)` must hold. dim, size(sample, 3 - dim), size(meang, 1) = "\
115 : //getStr([dim, size(sample, 3 - dim, IK), size(meang, 1, IK)]))
116 : #define CHECK_LEN_MEAN(PROC) \
117 : CHECK_ASSERTION(__LINE__, size(sample, 3 - dim, IK) == size(mean, 1, IK), \
118 : PROC//SK_": The condition `size(sample, 3 - dim) == size(mean)` must hold. dim, size(sample, 3 - dim), size(mean, 1) = "\
119 : //getStr([dim, size(sample, 3 - dim, IK), size(mean, 1, IK)]))
120 : #define CHECK_LEN_WEI(PROC,DIM) \
121 : CHECK_ASSERTION(__LINE__, size(sample, DIM, IK) == size(weight, 1, IK), \
122 : PROC//SK_": The condition `size(sample, dim) == size(weight)` must hold. dim, size(sample, dim), size(weight, 1) = "\
123 : //getStr([DIM, size(sample, DIM, IK), size(weight, 1, IK)]))
124 : #define CHECK_WEISUM(PROC) \
125 : CHECK_ASSERTION(__LINE__, abs(weisum - sum(weight)) < 1000 * epsilon(0._TKC), \
126 : PROC//SK_": The condition `0 < sum(weight)` must hold. weisum, sum(weight) = "//getStr([weisum, sum(weight)]))
127 :
128 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 : #if getCovMerged_ENABLED && New_ENABLED
130 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 :
132 : integer(IK) :: idim
133 460 : call setCovMerged(cov, covB, covA, meanDiff, fracA, uppDia)
134 : !do concurrent(idim = 2 : size(cov, 1, IK))
135 1314 : do idim = 2, size(cov, 1, IK)
136 2987 : cov(idim, 1 : idim - 1) = GET_CONJG(cov(1 : idim - 1, idim))
137 : end do
138 :
139 : !%%%%%%%%%%%%%%%%%%%
140 : #elif setCovMerged_ENABLED
141 : !%%%%%%%%%%%%%%%%%%%
142 :
143 : integer(IK) :: idim, jdim, ndim
144 : real(TKC) :: fracB, fracAB
145 2150 : fracB = 1._TKC - fracA
146 : ! Define the output value for setCovMerged.
147 : #if Old_ENABLED
148 : #define cov covB
149 : #elif New_ENABLED
150 14520 : CHECK_ASSERTION(__LINE__, all(shape(cov, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(shape(cov) == shape(covA))` must hold. shape(cov), shape(covA) = "//getStr([shape(cov, IK), shape(covA, IK)]))
151 : #else
152 : #error "Unrecognized interface."
153 : #endif
154 2150 : CHECK_ASSERTION(__LINE__, 0._TKC < fracA .and. fracA < 1._TKC, SK_"@setCovMerged(): The condition `0 < fracA .and. fracA < 1` must hold. fracA = "//getStr(fracA))
155 23650 : CHECK_ASSERTION(__LINE__, all(shape(covB, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(shape(covB, IK) == shape(covA, IK))` must hold. shape(covB), shape(covA) = "//getStr([shape(covB, IK), shape(covA, IK)]))
156 17200 : CHECK_ASSERTION(__LINE__, all(size(meanDiff, 1, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(size(meanDiff) == shape(covA))` must hold. size(meanDiff), shape(covA) = "//getStr([size(meanDiff, 1, IK), shape(covA, IK)]))
157 2150 : fracAB = fracA * fracB
158 2150 : ndim = size(cov, 1, IK)
159 8470 : do jdim = COL_RANGE(1_IK,ndim)
160 6320 : cov(jdim, jdim) = fracB * GET_RE(covB(jdim, jdim)) + fracA * GET_RE(covA(jdim, jdim)) + fracAB * GET_ABSQ(meanDiff(jdim))
161 16758 : do idim = ROW_RANGE(1_IK,jdim,ndim)
162 14608 : cov(idim, jdim) = fracB * covB(idim, jdim) + fracA * covA(idim, jdim) + fracAB * (meanDiff(idim) * GET_CONJG(meanDiff(jdim)))
163 : end do
164 : end do
165 : #undef cov
166 :
167 : !%%%%%%%%%%%%%%%%%%%%%%%
168 : #elif setCovMeanMerged_ENABLED
169 : !%%%%%%%%%%%%%%%%%%%%%%%
170 :
171 : TYPE_OF_SAMPLE :: idiff, temp
172 : integer(IK) :: idim, jdim, ndim
173 : real(TKC) :: fracB, fracAB
174 52502 : fracB = 1._TKC - fracA
175 : ! Define the output value for setCovMerged.
176 : #if Old_ENABLED
177 : #define mean meanB
178 : #define cov covB
179 : #elif New_ENABLED
180 568392 : CHECK_ASSERTION(__LINE__, all(shape(cov, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(shape(cov) == shape(covA))` must hold. shape(cov), shape(covA) = "//getStr([shape(cov, IK), shape(covA, IK)]))
181 413376 : CHECK_ASSERTION(__LINE__, all(size(mean, 1, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(size(mean) == shape(covA))` must hold. size(mean), shape(covA) = "//getStr([size(mean, 1, IK), shape(covA, IK)]))
182 : #else
183 : #error "Unrecognized interface."
184 : #endif
185 52502 : CHECK_ASSERTION(__LINE__, 0._TKC < fracA .and. fracA < 1._TKC, SK_"@setCovMerged(): The condition `0 < fracA .and. fracA < 1` must hold. fracA = "//getStr(fracA))
186 420016 : CHECK_ASSERTION(__LINE__, all(size(meanA, 1, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(size(meanA) == shape(covA))` must hold. size(meanA), shape(covA) = "//getStr([size(meanA, 1, IK), shape(covA, IK)]))
187 420016 : CHECK_ASSERTION(__LINE__, all(size(meanB, 1, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(size(meanB) == shape(covA))` must hold. size(meanB), shape(covA) = "//getStr([size(meanB, 1, IK), shape(covA, IK)]))
188 577522 : CHECK_ASSERTION(__LINE__, all(shape(covB, IK) == shape(covA, IK)), SK_"@setCovMerged(): The condition `all(shape(covB, IK) == shape(covA, IK))` must hold. shape(covB), shape(covA) = "//getStr([shape(covB, IK), shape(covA, IK)]))
189 52502 : fracAB = fracA * fracB
190 52502 : ndim = size(cov, 1, IK)
191 218025 : do jdim = COL_RANGE(1_IK,ndim)
192 165523 : temp = meanA(jdim) - meanB(jdim)
193 165523 : cov(jdim, jdim) = fracB * GET_RE(covB(jdim, jdim)) + fracA * GET_RE(covA(jdim, jdim)) + fracAB * GET_ABSQ(temp)
194 2484 : SET_CONJG(temp)
195 443330 : do idim = ROW_RANGE(1_IK,jdim,ndim)
196 225305 : idiff = meanA(idim) - meanB(idim)
197 390828 : cov(idim, jdim) = fracB * covB(idim, jdim) + fracA * covA(idim, jdim) + fracAB * idiff * temp
198 : end do
199 : end do
200 : !do concurrent(idim = 1 : size(covA, 1, IK))
201 218025 : do idim = 1, size(covA, 1, IK)
202 218025 : mean(idim) = fracA * meanA(idim) + fracB * meanB(idim)
203 : end do
204 : #undef mean
205 : #undef cov
206 :
207 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
208 : #elif getCov_ENABLED && CorStd_ENABLED
209 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210 :
211 : integer(IK) :: idim
212 1272 : call setCov(cov, uppDia, cor, subsetr, std)
213 : !do concurrent(idim = 2 : size(cov, 1, IK))
214 3820 : do idim = 2, size(cov, 1, IK)
215 8978 : cov(idim, 1 : idim - 1) = GET_CONJG(cov(1 : idim - 1, idim))
216 : end do
217 :
218 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
219 : #elif getCov_ENABLED && CorStd_ENABLED && 0
220 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 :
222 : #if ULD_ULD_ENABLED
223 : integer(IK) :: idim
224 : type(uppDia_type), parameter :: subset = uppDia_type(), subsetr = uppDia_type()
225 : #elif !(UXD_UXD_ENABLED || UXD_XLD_ENABLED || XLD_UXD_ENABLED || XLD_XLD_ENABLED)
226 : #error "Unrecognized interface."
227 : #endif
228 : call setCov(cov, subset, cor, subsetr, std)
229 : #if ULD_ULD_ENABLED
230 : !do concurrent(idim = 2 : size(cov, 1, IK))
231 : do idim = 2, size(cov, 1, IK)
232 : cov(idim, 1 : idim - 1) = GET_CONJG(cov(1 : idim - 1, idim))
233 : end do
234 : #endif
235 :
236 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
237 : #elif setCov_ENABLED && CorStd_ENABLED
238 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
239 :
240 : integer(IK) :: idim, jdim, ndim
241 6892 : ndim = size(cov, 1, IK)
242 6892 : CHECK_ASSERTION(__LINE__, ndim == size(cov, 2, IK), SK_"@setCov(): The condition `size(cov, 1) == size(cov, 2)` must hold. shape(cov) = "//getStr(shape(cov, IK)))
243 20676 : CHECK_ASSERTION(__LINE__, ndim == size(std, 1, IK), SK_"@setCov(): The condition `size(cov, 1) == size(std)` must hold. size(cov, 1), size(std) = "//getStr([ndim, size(std, 1, IK)]))
244 55136 : CHECK_ASSERTION(__LINE__, all(ndim == shape(cor, IK)), SK_"@setCov(): The condition `all(size(cov, 1) == shape(cor))` must hold. size(cov, 1), shape(cor) = "//getStr([ndim, shape(cor, IK)]))
245 27637 : CHECK_ASSERTION(__LINE__, all(0._TKC < std), SK_"@setCov(): The condition `all(0. < std)` must hold. std = "//getStr(std))
246 6892 : if (ndim == 0_IK) return
247 27635 : do jdim = COL_RANGE(1,ndim)
248 20745 : cov(jdim, jdim) = std(jdim)**2
249 55476 : do idim = ROW_RANGE(1,jdim,ndim)
250 : #if UXD_UXD_ENABLED || XLD_XLD_ENABLED
251 28435 : cov(idim, jdim) = std(idim) * std(jdim) * cor(idim, jdim)
252 : #elif UXD_XLD_ENABLED || XLD_UXD_ENABLED
253 20151 : cov(idim, jdim) = std(idim) * std(jdim) * GET_CONJG(cor(jdim, idim))
254 : #else
255 : #error "Unrecognized interface."
256 : #endif
257 : end do
258 : end do
259 : !#if UXD_UXD_ENABLED
260 : ! do jdim = 1_IK, ndim
261 : ! do idim = 1_IK, jdim - 1_IK
262 : ! cov(idim, jdim) = cor(idim, jdim) * std(idim) * std(jdim)
263 : ! end do
264 : ! cov(jdim, jdim) = std(jdim)**2
265 : ! end do
266 : !#elif UXD_XLD_ENABLED
267 : ! cov(1, 1) = std(1)**2
268 : ! do jdim = 2_IK, ndim
269 : ! do idim = 1_IK, jdim - 1_IK
270 : ! cov(idim, jdim) = GET_CONJG(cor(jdim, idim)) * std(idim) * std(jdim)
271 : ! end do
272 : ! cov(jdim, jdim) = std(jdim)**2
273 : ! end do
274 : !#elif XLD_UXD_ENABLED
275 : ! cov(1, 1) = std(1)**2
276 : ! do jdim = 2_IK, ndim
277 : ! do idim = 1_IK, jdim - 1_IK
278 : ! cov(jdim, idim) = GET_CONJG(cor(idim, jdim)) * std(idim) * std(jdim)
279 : ! end do
280 : ! cov(jdim, jdim) = std(jdim)**2
281 : ! end do
282 : !#elif XLD_XLD_ENABLED
283 : ! cov(1, 1) = std(1)**2
284 : ! do jdim = 2_IK, ndim
285 : ! do idim = 1_IK, jdim - 1_IK
286 : ! cov(jdim, idim) = cor(jdim, idim) * std(idim) * std(jdim)
287 : ! end do
288 : ! cov(jdim, jdim) = std(jdim)**2
289 : ! end do
290 : !#else
291 : !#error "Unrecognized interface."
292 : !#endif
293 :
294 : !%%%%%%%%%%%%%
295 : #elif getCov_ENABLED
296 : !%%%%%%%%%%%%%
297 :
298 : type(uppDia_type), parameter :: subset = uppDia_type()
299 : integer(IK) :: ndim, idim, nsam
300 : real(TKC) :: normfac
301 : #if WNO_ENABLED
302 : #define WEIGHT_ARGS
303 : #elif WTI_ENABLED || WTR_ENABLED
304 : #define WEIGHT_ARGS, weight, weisum
305 : TYPE_OF_WEIGHT :: weisum
306 : #else
307 : #error "Unrecognized interface."
308 : #endif
309 : #if XY_ENABLED
310 : TYPE_OF_SAMPLE :: mean(2)
311 7218 : call setCovMean(cov, mean, x, y WEIGHT_ARGS, [x(1), y(1)])
312 802 : nsam = size(x, 1, IK)
313 : #elif ULD_ENABLED
314 12564 : TYPE_OF_SAMPLE, dimension(size(sample, 3 - dim, IK)) :: mean, meang
315 4494 : nsam = size(sample, dim, IK)
316 6282 : if (dim == 1_IK) then
317 5058 : meang = sample(1,:)
318 : else
319 19803 : meang = sample(:,1)
320 : end if
321 6282 : call setCovMean(cov, subset, mean, sample, dim WEIGHT_ARGS, meang)
322 : #else
323 : #error "Unrecognized interface."
324 : #endif
325 : ! Symmetrize.
326 6282 : ndim = size(cov, 1, IK)
327 : !do concurrent(idim = 1 : ndim)
328 32079 : do idim = 1, ndim
329 58968 : cov(idim, 1 : idim - 1) = GET_CONJG(cov(1 : idim - 1, idim))
330 : end do
331 : ! Correct if requested.
332 8688 : if (present(correction)) then
333 2400 : CHECK_ASSERTION(__LINE__, same_type_as(correction, fweight) .or. same_type_as(correction, rweight), SK_"@getCov(): The condition `same_type_as(correction, fweight) .or. same_type_as(correction, rweight)` must hold.")
334 : #if WNO_ENABLED
335 800 : normfac = getVarCorrection(real(nsam, TKC))
336 : #elif (WTI_ENABLED || WTR_ENABLED)
337 1600 : if (same_type_as(correction, fweight)) then
338 800 : normfac = getVarCorrection(real(weisum, TKC))
339 800 : elseif (same_type_as(correction, rweight)) then
340 4448 : normfac = getVarCorrection(real(weisum, TKC), real(sum(weight**2), TKC))
341 : end if
342 : #else
343 : #error "Unrecognized interface."
344 : #endif
345 25896 : cov = cov * normfac
346 : end if
347 : #undef WEIGHT_ARGS
348 : #undef GET_SAM
349 :
350 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
351 : #elif setCov_ENABLED && XY_ENABLED
352 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
353 :
354 : integer(IK) :: isam, nsam
355 : #if WTI_ENABLED || WTR_ENABLED
356 : TYPE_OF_SAMPLE :: temp
357 : #endif
358 : TYPE_OF_SAMPLE :: cxy
359 : real(TKC) :: cxx, cyy
360 : real(TKC) :: normFac
361 : #if Avg_ENABLED
362 : TYPE_OF_SAMPLE :: difx, dify
363 13218 : CHECK_ASSERTION(__LINE__, size(mean, 1, IK) == 2_IK, SK_"@setCov(): The condition `size(mean) == 2` must hold. size(mean) = "//getStr(size(mean, 1, IK)))
364 : #elif !Org_ENABLED
365 : #error "Unrecognized interface."
366 : #endif
367 15624 : nsam = size(x, 1, IK)
368 46872 : CHECK_ASSERTION(__LINE__, size(x, 1, IK) == size(y, 1, IK), SK_"@setCov(): The condition `size(x) == size(y)` must hold. size(x), size(y) = "//getStr([size(x, 1, IK), size(y, 1, IK)]))
369 46872 : CHECK_ASSERTION(__LINE__, all(shape(cov, IK) == 2_IK), SK_"@setCov(): The condition `all(shape(cov) == 2)` must hold. shape(cov) = "//getStr(shape(cov, IK)))
370 15624 : CHECK_ASSERTION(__LINE__, 1_IK < nsam, SK_"@setCov(): The condition `1 < size(x)` must hold. size(x) = "//getStr(nsam))
371 31490 : CHECK_ASSERTION(__LINE__, any(x(1) /= x), SK_"@setCov(): The condition `any(x(1) /= x)` must hold. x = "//getStr(x))
372 31477 : CHECK_ASSERTION(__LINE__, any(y(1) /= y), SK_"@setCov(): The condition `any(y(1) /= y)` must hold. y = "//getStr(y))
373 : #if WNO_ENABLED
374 3460 : normFac = 1._TKC / real(nsam, TKC)
375 : #elif WTI_ENABLED || WTR_ENABLED
376 31242 : CHECK_ASSERTION(__LINE__, size(x, 1, IK) == size(weight, 1, IK), SK_"@setCov(): The condition `size(x) == size(weight)` must hold. size(x), size(weight) = "//getStr([size(x, 1, IK), size(weight, 1, IK)]))
377 301122 : CHECK_ASSERTION(__LINE__, all(0 <= weight), SK_"@setCov(): The condition `all(0 <= weight)` must hold. pack(weight, weight < 0) = "//getStr(pack(weight, weight < 0)))
378 6914 : normFac = 1._TKC / real(weisum, TKC)
379 : #else
380 : #error "Unrecognized interface."
381 : #endif
382 4200 : cxx = 0._TKC
383 1050 : cxy = 0._TKC
384 4200 : cyy = 0._TKC
385 233685 : do isam = 1, nsam
386 : #if Avg_ENABLED && (WTI_ENABLED || WTR_ENABLED)
387 138048 : difx = x(isam) - mean(1)
388 138048 : dify = y(isam) - mean(2)
389 138048 : temp = dify * weight(isam)
390 138048 : cyy = cyy + GET_PROD(dify,temp)
391 138048 : cxy = cxy + difx * GET_CONJG(temp)
392 146858 : cxx = cxx + GET_ABSQ(difx) * weight(isam)
393 : #elif Org_ENABLED && (WTI_ENABLED || WTR_ENABLED)
394 7306 : temp = y(isam) * weight(isam)
395 7306 : cyy = cyy + GET_PROD(y(isam),temp)
396 7306 : cxy = cxy + x(isam) * GET_CONJG(temp)
397 8910 : cxx = cxx + GET_ABSQ(x(isam)) * weight(isam)
398 : #elif Avg_ENABLED && WNO_ENABLED
399 69054 : difx = x(isam) - mean(1)
400 69054 : dify = y(isam) - mean(2)
401 69054 : cyy = cyy + GET_ABSQ(dify)
402 69054 : cxy = cxy + difx * GET_CONJG(dify)
403 73462 : cxx = cxx + GET_ABSQ(difx)
404 : #elif Org_ENABLED && WNO_ENABLED
405 3653 : cyy = cyy + GET_ABSQ(y(isam))
406 3653 : cxy = cxy + x(isam) * GET_CONJG(y(isam))
407 4455 : cxx = cxx + GET_ABSQ(x(isam))
408 : #else
409 : #error "Unrecognized interface."
410 : #endif
411 : end do
412 15624 : cov(1,1) = normFac * cxx
413 15624 : cov(2,2) = normFac * cyy
414 15624 : cov(1,2) = normFac * cxy
415 15624 : cov(2,1) = GET_CONJG(cov(1,2))
416 :
417 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
418 : #elif setCov_ENABLED && (UXD_ENABLED || XLD_ENABLED)
419 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
420 :
421 : real(TKC) :: normFac
422 : TYPE_OF_SAMPLE :: temp
423 : integer(IK) :: idim, jdim, isam, ndim, nsam
424 : #if Avg_ENABLED
425 : #if WTI_ENABLED || WTR_ENABLED
426 : TYPE_OF_SAMPLE :: diff
427 : #elif !WNO_ENABLED
428 : #error "Unrecognized interface."
429 : #endif
430 281608 : CHECK_LEN_MEAN(SK_"@setCov()")
431 : #elif !Org_ENABLED
432 : #error "Unrecognized interface."
433 : #endif
434 78828 : ndim = size(sample, 3 - dim, IK)
435 78828 : nsam = size(sample, dim, IK)
436 236484 : CHECK_VAL_DIM(SK_"@setCov()")
437 709452 : CHECK_SHAPE_COV(SK_"@setCov()")
438 472968 : CHECK_VAL_NSAM(SK_"@setCov()",dim)
439 472968 : CHECK_VAL_NDIM(SK_"@setCov()",3 - dim)
440 : #if WNO_ENABLED
441 9330 : normFac = 1._TKC / real(nsam, TKC)
442 : #elif WTI_ENABLED || WTR_ENABLED
443 2058380 : CHECK_WEISUM(SK_"@setCov()")
444 994441 : CHECK_SUM_WEI(SK_"@setCov()")
445 994441 : CHECK_VAL_WEI(SK_"@setCov()")
446 277992 : CHECK_LEN_WEI(SK_"@setCov()",dim)
447 69498 : normFac = 1._TKC / real(weisum, TKC)
448 : #else
449 : #error "Unrecognized interface."
450 : #endif
451 : ! Compute the cov.
452 78828 : if (dim == 2_IK) then
453 264835 : do jdim = COL_RANGE(1_IK,ndim)
454 803153 : CHECK_ASSERTION(__LINE__, any(sample(jdim,1) /= sample(jdim,:)), SK_"@setCov(): The condition `any(sample(jdim,1) /= sample(jdim,:))` must hold. jdim = "//getStr(jdim)//SK_", sample(jdim,:) = "//getStr(reshape(sample(jdim,:),[size(sample,2),1])))
455 200668 : cov(jdim, jdim) = 0._TKC
456 537429 : do idim = ROW_RANGE(1_IK,jdim,ndim)
457 473262 : cov(idim, jdim) = 0._TKC
458 : end do
459 : end do
460 931187 : do isam = 1, nsam
461 3973987 : do jdim = COL_RANGE(1_IK,ndim)
462 : #if Avg_ENABLED && (WTI_ENABLED || WTR_ENABLED)
463 2846724 : diff = sample(jdim, isam) - mean(jdim)
464 2846724 : temp = diff * weight(isam)
465 2846724 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_PROD(diff,temp)
466 : #elif Org_ENABLED && (WTI_ENABLED || WTR_ENABLED)
467 32588 : temp = sample(jdim, isam) * weight(isam)
468 32588 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_PROD(sample(jdim, isam),temp)
469 : #elif Avg_ENABLED && WNO_ENABLED
470 147064 : temp = sample(jdim, isam) - mean(jdim)
471 147064 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_PROD(temp,temp)
472 : #elif Org_ENABLED && WNO_ENABLED
473 16424 : temp = sample(jdim, isam)
474 16424 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_PROD(sample(jdim, isam),temp)
475 : #else
476 : #error "Unrecognized interface."
477 : #endif
478 59361 : SET_CONJG(temp)
479 8352313 : do idim = ROW_RANGE(1_IK,jdim,ndim)
480 7485293 : cov(idim, jdim) = cov(idim, jdim) + GET_SHIFTED(sample(idim, isam),mean(idim)) * temp
481 : end do
482 : end do
483 : end do
484 264835 : do jdim = COL_RANGE(1_IK,ndim)
485 200668 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) * normFac
486 537429 : do idim = ROW_RANGE(1_IK,jdim,ndim)
487 473262 : cov(idim, jdim) = cov(idim, jdim) * normFac
488 : end do
489 : end do
490 : else ! dim = 1_IK
491 57219 : do jdim = COL_RANGE(1_IK,ndim)
492 85479 : CHECK_ASSERTION(__LINE__, any(sample(1,jdim) /= sample(:,jdim)), SK_"@setCov(): The condition `any(sample(1,jdim) /= sample(:,jdim))` must hold. jdim = "//getStr(jdim)//SK_", sample(:,jdim) = "//getStr(sample(:,jdim)))
493 42558 : cov(jdim, jdim) = 0._TKC
494 545038 : do isam = 1, nsam
495 : #if Avg_ENABLED
496 439849 : temp = sample(isam, jdim) - mean(jdim)
497 468989 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_WEIGHTED(GET_ABSQ(temp),weight(isam))
498 : #elif Org_ENABLED
499 76049 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_WEIGHTED(GET_ABSQ(sample(isam, jdim)),weight(isam))
500 : #endif
501 : end do
502 42558 : cov(jdim, jdim) = cov(jdim, jdim) * normFac
503 112086 : do idim = ROW_RANGE(1_IK,jdim,ndim)
504 54867 : cov(idim, jdim) = 0._TKC
505 723626 : do isam = 1, nsam
506 723626 : cov(idim, jdim) = cov(idim, jdim) + GET_WEIGHTED(GET_SHIFTED(sample(isam, idim),mean(idim)) * GET_CONJG(GET_SHIFTED(sample(isam, jdim),mean(jdim))),weight(isam))
507 : end do
508 97425 : cov(idim, jdim) = cov(idim, jdim) * normFac
509 : end do
510 : end do
511 : end if ! dim
512 :
513 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
514 : #elif setCovMean_ENABLED && XY_ENABLED
515 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
516 :
517 : integer(IK) :: isam
518 : TYPE_OF_SAMPLE :: cxy, difx, dify
519 : real(TKC) :: cxx, cyy
520 : real(TKC) :: normFac
521 : #if WNO_ENABLED
522 : real(TKC) :: weisum
523 1603 : weisum = size(x, 1, IK)
524 : #elif WTI_ENABLED || WTR_ENABLED
525 : TYPE_OF_SAMPLE :: temp
526 9618 : CHECK_ASSERTION(__LINE__, size(x, 1, IK) == size(weight, 1, IK), SK_"@setCovMean(): The condition `size(x) == size(weight)` must hold. size(x), size(weight) = "//getStr([size(x, 1, IK), size(weight, 1, IK)]))
527 32514 : CHECK_ASSERTION(__LINE__, all(0 <= weight), SK_"@setCovMean(): The condition `all(0 <= weight)` must hold. pack(weight, weight < 0) = "//getStr(pack(weight, weight < 0)))
528 3206 : weisum = ZEROW
529 : #else
530 : #error "Unrecognized interface."
531 : #endif
532 4809 : CHECK_ASSERTION(__LINE__, size(mean, 1, IK) == 2_IK, SK_"@setCovMean(): The condition `size(mean) == 2` must hold. size(mean) = "//getStr(size(mean, 1, IK)))
533 4809 : CHECK_ASSERTION(__LINE__, size(meang, 1, IK) == 2_IK, SK_"@setCovMean(): The condition `size(meang) == 2` must hold. size(meang) = "//getStr(size(meang, 1, IK)))
534 14427 : CHECK_ASSERTION(__LINE__, all(shape(cov, IK) == 2_IK), SK_"@setCovMean(): The condition `all(shape(cov) == 2)` must hold. shape(cov) = "//getStr(shape(cov, IK)))
535 14427 : CHECK_ASSERTION(__LINE__, size(x, 1, IK) == size(y, 1, IK), SK_"@setCovMean(): The condition `size(x) == size(y)` must hold. size(x), size(y) = "//getStr([size(x, 1, IK), size(y, 1, IK)]))
536 94977 : CHECK_ASSERTION(__LINE__, all([minval(x) <= meang(1) .and. meang(2) <= maxval(y)]), SK_"@setCovMean(): The condition `all([minval(x) <= meang(1) .and. meang(2) <= maxval(y)])` must hold. minval(x), meang, maxval(y) = "//getStr([minval(x), meang, maxval(y)]))
537 4809 : CHECK_ASSERTION(__LINE__, 1_IK < size(x, 1, IK), SK_"@setCovMean(): The condition `1 < size(x)` must hold. size(x) = "//getStr(size(x, 1, IK)))
538 9618 : CHECK_ASSERTION(__LINE__, any(x(1) /= x), SK_"@setCovMean(): The condition `any(x(1) /= x)` must hold. x = "//getStr(x))
539 9622 : CHECK_ASSERTION(__LINE__, any(y(1) /= y), SK_"@setCovMean(): The condition `any(y(1) /= y)` must hold. y = "//getStr(y))
540 2400 : cxx = 0._TKC
541 2400 : cyy = 0._TKC
542 3000 : cxy = 0._TKC
543 14427 : mean = 0._TKC
544 26790 : do isam = 1, size(x, 1, IK)
545 21981 : difx = x(isam) - meang(1)
546 21981 : dify = y(isam) - meang(2)
547 : #if WTI_ENABLED || WTR_ENABLED
548 14654 : weisum = weisum + weight(isam)
549 14654 : temp = difx * weight(isam)
550 14654 : mean(1) = mean(1) + temp
551 14654 : cxx = cxx + GET_PROD(difx,temp)
552 14654 : temp = dify * weight(isam)
553 14654 : mean(2) = mean(2) + temp
554 14654 : cyy = cyy + GET_PROD(dify,temp)
555 17860 : cxy = cxy + difx * GET_CONJG(temp)
556 : #elif WNO_ENABLED
557 7327 : mean(1) = mean(1) + difx
558 7327 : cxx = cxx + GET_ABSQ(difx)
559 7327 : mean(2) = mean(2) + dify
560 7327 : cyy = cyy + GET_ABSQ(dify)
561 8930 : cxy = cxy + difx * GET_CONJG(dify)
562 : #else
563 : #error "Unrecognized interface."
564 : #endif
565 : end do
566 2406 : normFac = 1._TKC / weisum
567 14427 : mean = mean * normFac
568 4809 : cov(1,1) = (cxx - GET_ABSQ(mean(1)) * weisum) * normFac
569 4809 : cov(2,2) = (cyy - GET_ABSQ(mean(2)) * weisum) * normFac
570 4809 : cov(1,2) = (cxy - mean(1) * GET_CONJG(mean(2)) * weisum) * normFac
571 4809 : cov(2,1) = GET_CONJG(cov(1,2))
572 14427 : mean = mean + meang
573 :
574 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
575 : #elif setCovMean_ENABLED && (UXD_ENABLED || XLD_ENABLED)
576 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
577 :
578 : real(TKC) :: normFac
579 : TYPE_OF_SAMPLE :: diff
580 : integer(IK) :: idim, jdim, isam, ndim, nsam
581 : #if WNO_ENABLED
582 : #define SET(X)
583 : #define temp diff
584 : #define INCREMENT(X,Y)
585 : real(TKC) :: weisum
586 5698 : weisum = size(sample, dim, IK)
587 : #elif WTI_ENABLED || WTR_ENABLED
588 : #define INCREMENT(X,Y)X = X + Y
589 : #define SET(X)X
590 : TYPE_OF_SAMPLE :: temp
591 16784 : CHECK_LEN_WEI(SK_"@setCovMean()",dim)
592 23652 : CHECK_VAL_WEI(SK_"@setCovMean()")
593 4196 : weisum = ZEROW
594 : #else
595 : #error "Unrecognized interface."
596 : #endif
597 9894 : nsam = size(sample, dim, IK)
598 9894 : ndim = size(sample, 3 - dim, IK)
599 39576 : CHECK_LEN_MEAN(SK_"@setCovMean()")
600 59364 : CHECK_VAL_NSAM(SK_"@setCovMean()",dim)
601 39054 : CHECK_VAL_MEANG(SK_"@setCovMean()",dim)
602 39576 : CHECK_LEN_MEANG(SK_"@setCovMean()")
603 89046 : CHECK_SHAPE_COV(SK_"@setCovMean()")
604 29682 : CHECK_VAL_DIM(SK_"@setCovMean()")
605 59364 : CHECK_VAL_NSAM(SK_"@setCovMean()",dim)
606 59364 : CHECK_VAL_NDIM(SK_"@setCovMean()",3 - dim)
607 : ! Compute the cov.
608 9894 : if (dim == 2_IK) then
609 26934 : do jdim = COL_RANGE(1_IK,ndim)
610 134884 : CHECK_ASSERTION(__LINE__, any(sample(jdim,1) /= sample(jdim,:)), SK_"@setCovMean(): The condition `any(sample(jdim,1) /= sample(jdim,:))` must hold. jdim = "//getStr(jdim)//SK_", sample(jdim,:) = "//getStr(sample(jdim,:)))
611 20145 : mean(jdim) = 0._TKC
612 20145 : cov(jdim, jdim) = 0._TKC
613 53520 : do idim = ROW_RANGE(1_IK,jdim,ndim)
614 46731 : cov(idim, jdim) = 0._TKC
615 : end do
616 : end do
617 38649 : do isam = 1, nsam
618 10004 : INCREMENT(weisum,weight(isam))
619 133239 : do jdim = COL_RANGE(1_IK,ndim)
620 94590 : diff = sample(jdim, isam) - meang(jdim)
621 29118 : SET(temp = diff * weight(isam))
622 94590 : mean(jdim) = mean(jdim) + temp
623 94590 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_PROD(diff,temp)
624 45163 : SET_CONJG(temp)
625 251485 : do idim = ROW_RANGE(1_IK,jdim,ndim)
626 219625 : cov(idim, jdim) = cov(idim, jdim) + (sample(idim, isam) - meang(idim)) * temp
627 : end do
628 : end do
629 : end do
630 3570 : normFac = 1._TKC / weisum
631 26934 : mean = mean * normFac
632 26934 : do jdim = COL_RANGE(1_IK,ndim)
633 20145 : GET_RE(cov(jdim, jdim)) = (GET_RE(cov(jdim, jdim)) - GET_ABSQ(mean(jdim)) * weisum) * normFac
634 14613 : temp = GET_CONJG(mean(jdim))
635 53520 : do idim = ROW_RANGE(1_IK,jdim,ndim)
636 46731 : cov(idim, jdim) = (cov(idim, jdim) - mean(idim) * temp * weisum) * normFac
637 : end do
638 : end do
639 26934 : mean = mean + meang
640 : else ! dim == 1
641 : ! Compute `weisum` and the first element of `mean` in the first round.
642 6214 : CHECK_ASSERTION(__LINE__, any(sample(1,FIRST) /= sample(:,FIRST)), SK_"@setCovMean(): The condition `any(sample(1,idim) /= sample(:,idim))` must hold. idim = "//getStr(FIRST)//SK_", sample(:,idim) = "//getStr(sample(:,FIRST)))
643 3105 : cov(FIRST, FIRST) = 0._TKC
644 3105 : mean(FIRST) = 0._TKC
645 17283 : do isam = 1, nsam
646 9452 : INCREMENT(weisum,weight(isam))
647 14178 : diff = sample(isam, FIRST) - meang(FIRST)
648 9452 : SET(temp = diff * weight(isam))
649 14178 : GET_RE(cov(FIRST, FIRST)) = GET_RE(cov(FIRST, FIRST)) + GET_PROD(diff,temp)
650 17283 : mean(FIRST) = mean(FIRST) + temp
651 : end do
652 3105 : normFac = 1._TKC / weisum
653 3105 : mean(FIRST) = mean(FIRST) * normFac
654 3105 : GET_RE(cov(FIRST, FIRST)) = (GET_RE(cov(FIRST, FIRST)) - GET_ABSQ(mean(FIRST)) * weisum) * normFac
655 : ! Compute the rest of the `mean` elements in the second round.
656 9015 : do idim = OFF_RANGE(1_IK,ndim,1_IK)
657 11820 : CHECK_ASSERTION(__LINE__, any(sample(1,idim) /= sample(:,idim)), SK_"@setCovMean(): The condition `any(sample(1,idim) /= sample(:,idim))` must hold. idim = "//getStr(idim)//SK_", sample(:,idim) = "//getStr(sample(:,idim)))
658 5910 : cov(idim, FIRST) = 0._TKC
659 5910 : mean(idim) = 0._TKC
660 32784 : do isam = 1, nsam
661 26874 : diff = sample(isam, idim) - meang(idim)
662 17916 : SET(temp = diff * weight(isam))
663 26874 : cov(idim, FIRST) = cov(idim, FIRST) + temp * GET_CONJG((sample(isam, FIRST) - meang(FIRST)))
664 32784 : mean(idim) = mean(idim) + temp
665 : end do
666 5910 : mean(idim) = mean(idim) * normFac
667 5910 : cov(idim, FIRST) = (cov(idim, FIRST) - mean(idim) * GET_CONJG(mean(FIRST)) * weisum) * normFac
668 9015 : mean(idim) = mean(idim) + meang(idim)
669 : end do
670 3105 : mean(FIRST) = mean(FIRST) + meang(FIRST) ! This normalization must be done right here and not any sooner.
671 : ! Now use the computed mean to calculate the rest of the covariance matrix using the normal algorithm.
672 9015 : do jdim = OFF_RANGE(1_IK,ndim,1_IK)
673 5910 : cov(jdim, jdim) = 0._TKC
674 32784 : do isam = 1, nsam
675 26874 : diff = sample(isam, jdim) - mean(jdim)
676 32784 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) + GET_WEIGHTED(GET_ABSQ(diff),weight(isam))
677 : end do
678 5910 : GET_RE(cov(jdim, jdim)) = GET_RE(cov(jdim, jdim)) * normFac
679 14829 : do idim = ROW_RANGE(1_IK,jdim,ndim)
680 5814 : cov(idim, jdim) = 0._TKC
681 32085 : do isam = 1, nsam
682 32085 : cov(idim, jdim) = cov(idim, jdim) + GET_WEIGHTED((sample(isam, idim) - mean(idim)) * GET_CONJG((sample(isam, jdim) - mean(jdim))),weight(isam))
683 : end do
684 11724 : cov(idim, jdim) = cov(idim, jdim) * normFac
685 : end do
686 : end do
687 : end if
688 : #undef INCREMENT
689 : #undef temp
690 : #undef SET
691 :
692 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
693 : #elif setCov_ENABLED && Wei_ENABLED && Prob_ENABLED
694 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
695 :
696 : #error "Unrecognized interface. This is a relic of the past from the era of frequentist ignorance."
697 : ! ! Define the looping ranges for the corresponding matrix subsets.
698 : !if (setCov_ENABLED || setCovMerged_ENABLED) && UXD_ENABLED
699 : !define COL_RANGE 2_IK, ndim
700 : !define ROW_RANGE 1_IK, jdim
701 : !define FIRST 1_IK
702 : !elif (setCov_ENABLED || setCovMerged_ENABLED) && XLD_ENABLED
703 : !define COL_RANGE ndim - 1_IK, 1_IK, -1_IK
704 : !define ROW_RANGE jdim, ndim
705 : !define FIRST ndim
706 : !elif setCovMerged_ENABLED || (setCov_ENABLED && !(XY_ENABLED || CorStd_ENABLED))
707 : !error "Unrecognized interface."
708 : !endif
709 : ! integer(IK) :: zeroWeightCount
710 : ! integer(IK) :: idim, jdim, isam,ndim,nsam
711 : ! real(TKC) :: normFac
712 : !#if Def_ENABLED
713 : ! TYPE_OF_WEIGHT :: weisum
714 : ! real(TKC), allocatable :: mean(:)
715 : !#elif Avg_ENABLED
716 : ! CHECK_LEN_MEAN(SK_"@setCov()")
717 : !#else
718 : !#error "Unrecognized interface."
719 : !#endif
720 : ! CHECK_VAL_DIM(SK_"@setCov()")
721 : ! CHECK_SHAPE_COV(SK_"@setCov()")
722 : ! CHECK_SUM_WEI(SK_"@setCov()")
723 : ! CHECK_VAL_WEI(SK_"@setCov()")
724 : ! CHECK_LEN_WEI(SK_"@setCov()",dim)
725 : ! CHECK_VAL_NSAM(SK_"@setCov()",dim)
726 : ! CHECK_VAL_NDIM(SK_"@setCov()"),3 - dim)
727 : ! nsam = size(sample, dim, IK)
728 : ! ndim = size(sample, 3 - dim, IK)
729 : ! zeroWeightCount = 0_IK
730 : ! ! Compute the cov.
731 : ! if (dim == 2_IK) then
732 : ! do jdim = 1, ndim; do idim = ROW_RANGE; cov(idim, jdim) = 0._TKC; end do; end do;
733 : !#if Def_ENABLED
734 : ! if (shifted) then
735 : ! ! probability-weight correction and shifted.
736 : ! weisum = ZEROW
737 : ! do isam = 1, nsam
738 : ! weisum = weisum + weight(isam)
739 : ! if (weight(isam) == 0._TKC) zeroWeightCount = zeroWeightCount + 1
740 : ! do jdim = 1, ndim
741 : ! do idim = ROW_RANGE
742 : ! cov(idim, jdim) = cov(idim, jdim) + sample(idim, isam) * sample(jdim, isam) * weight(isam)**2
743 : ! end do
744 : ! end do
745 : ! end do
746 : ! normFac = real(nsam - zeroWeightCount, TKC) / (real(nsam - zeroWeightCount - 1, TKC) * weisum**2)
747 : ! do jdim = 1, ndim; do idim = ROW_RANGE; cov(idim, jdim) = cov(idim, jdim) * normFac; end do; end do;
748 : ! return
749 : ! end if
750 : ! allocate(mean(ndim))
751 : ! call setMean(mean, sample, dim, weight, weisum)
752 : !#endif
753 : ! ! probability-weight correction and not shifted.
754 : ! do isam = 1, nsam
755 : ! if (weight(isam) == 0._TKC) zeroWeightCount = zeroWeightCount + 1
756 : ! do jdim = 1, ndim
757 : ! do idim = ROW_RANGE
758 : ! cov(idim, jdim) = cov(idim, jdim) + (sample(idim, isam) - mean(idim)) * (sample(jdim, isam) - mean(jdim)) * weight(isam)**2
759 : ! end do
760 : ! end do
761 : ! end do
762 : ! normFac = real(nsam - zeroWeightCount, TKC) / (real(nsam - zeroWeightCount - 1, TKC) * weisum**2)
763 : ! do jdim = 1, ndim; do idim = ROW_RANGE; cov(idim, jdim) = cov(idim, jdim) * normFac; end do; end do;
764 : ! else ! dim = 1_IK
765 : !#if Def_ENABLED
766 : ! if (shifted) then
767 : ! ! probability-weight correction and shifted.
768 : ! weisum = ZEROW
769 : ! cov(FIRST, FIRST) = 0._TKC
770 : ! do isam = 1, nsam
771 : ! weisum = weisum + weight(isam)
772 : ! cov(FIRST, FIRST) = cov(FIRST, FIRST) + (sample(isam, FIRST) * weight(isam))**2
773 : ! if (weight(isam) == 0._TKC) zeroWeightCount = zeroWeightCount + 1
774 : ! end do
775 : ! normFac = real(nsam - zeroWeightCount, TKC) / (real(nsam - zeroWeightCount - 1, TKC) * weisum**2)
776 : ! cov(FIRST, FIRST) = cov(FIRST, FIRST) * normFac
777 : ! do jdim = COL_RANGE
778 : ! do idim = ROW_RANGE
779 : ! cov(idim, jdim) = 0._TKC
780 : ! do isam = 1, nsam
781 : ! cov(idim, jdim) = cov(idim, jdim) + sample(isam, idim) * sample(isam, jdim) * weight(isam)**2
782 : ! end do
783 : ! cov(idim, jdim) = cov(idim, jdim) * normFac
784 : ! end do
785 : ! end do
786 : ! return
787 : ! end if
788 : ! allocate(mean(ndim))
789 : ! call setMean(mean, sample, dim, weight, weisum)
790 : !#endif
791 : ! ! probability-weight correction and not shifted.
792 : ! cov(FIRST, FIRST) = 0._TKC
793 : ! do isam = 1, nsam
794 : ! if (weight(isam) == 0._TKC) zeroWeightCount = zeroWeightCount + 1
795 : ! cov(FIRST, FIRST) = cov(FIRST, FIRST) + ((sample(isam, FIRST) - mean(FIRST)) * weight(isam))**2
796 : ! end do
797 : ! normFac = real(nsam - zeroWeightCount, TKC) / (real(nsam - zeroWeightCount - 1, TKC) * weisum**2)
798 : ! cov(FIRST, FIRST) = cov(FIRST, FIRST) * normFac
799 : ! do jdim = COL_RANGE
800 : ! do idim = ROW_RANGE
801 : ! cov(idim, jdim) = 0._TKC
802 : ! do isam = 1, nsam
803 : ! cov(idim, jdim) = cov(idim, jdim) + (sample(isam, idim) - mean(idim) * (sample(isam, jdim) - mean(jdim) * weight(isam)**2
804 : ! end do
805 : ! cov(idim, jdim) = cov(idim, jdim) * normFac
806 : ! end do
807 : ! end do
808 : ! end if ! dim
809 : !#if Def_ENABLED
810 : ! deallocate(mean) ! Bypass gfortran bug for automatic deallocation of local heap arrays.
811 : !#endif
812 : !undef COL_RANGE
813 :
814 : #else
815 : !%%%%%%%%%%%%%%%%%%%%%%%%
816 : #error "Unrecognized interface."
817 : !%%%%%%%%%%%%%%%%%%%%%%%%
818 : #endif
819 : #undef TYPE_OF_WEIGHT
820 : #undef TYPE_OF_SAMPLE
821 : #undef CHECK_SHAPE_COV
822 : #undef CHECK_LEN_MEANG
823 : #undef CHECK_VAL_MEANG
824 : #undef CHECK_LEN_MEAN
825 : #undef CHECK_VAL_NSAM
826 : #undef CHECK_VAL_NDIM
827 : #undef CHECK_LEN_WEI
828 : #undef CHECK_VAL_DIM
829 : #undef CHECK_SUM_WEI
830 : #undef CHECK_VAL_WEI
831 : #undef CHECK_WEISUM
832 : #undef GET_WEIGHTED
833 : #undef GET_SHIFTED
834 : #undef GET_CONJG
835 : #undef SET_CONJG
836 : #undef OFF_RANGE
837 : #undef COL_RANGE
838 : #undef ROW_RANGE
839 : #undef GET_ABSQ
840 : #undef GET_PROD
841 : #undef GET_RE
842 : #undef FIRST
843 : #undef ZEROW
|
