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 : ! Set the conjugation rule.
28 : #if CK_ENABLED
29 : #define TYPE_OF_SAMPLE complex(TKC)
30 : #define GET_ABSQ(X)(real(X)**2 + aimag(X)**2)
31 : #define GET_PROD(X,Y)(X%re * Y%re + X%im * Y%im)
32 : #elif RK_ENABLED
33 : #define TYPE_OF_SAMPLE real(TKC)
34 : #define GET_PROD(X,Y)X * Y
35 : #define GET_ABSQ(X)X**2
36 : #else
37 : #error "Unrecognized interface."
38 : #endif
39 : ! Set the shifting rule.
40 : #if setVar_ENABLED && Org_ENABLED
41 : #define GET_SHIFTED(X,Y)X
42 : #elif setVar_ENABLED && Avg_ENABLED
43 : #define GET_SHIFTED(X,Y)(X - Y)
44 : #elif setVar_ENABLED
45 : #error "Unrecognized interface."
46 : #endif
47 : ! Set the weighting rule.
48 : #if setVar_ENABLED && WNO_ENABLED
49 : #define GET_WEIGHTED(X,W)X
50 : #elif setVar_ENABLED && (WTI_ENABLED || WTR_ENABLED)
51 : #define GET_WEIGHTED(X,W)X * W
52 : #elif setVar_ENABLED
53 : #error "Unrecognized interface."
54 : #endif
55 : ! Define weight type and kind.
56 : #if WTI_ENABLED
57 : #define TYPE_OF_WEIGHT integer(IK)
58 : #elif WTR_ENABLED
59 : #define TYPE_OF_WEIGHT real(TKC)
60 : #elif !WNO_ENABLED && (getVar_ENABLED || setVarMean_ENABLED)
61 : #error "Unrecognized interface."
62 : #endif
63 : ! Define runtime sanity checks.
64 : #define CHECK_VAL_NSAM(PROC,DIM) \
65 : CHECK_ASSERTION(__LINE__, 1 < size(sample, DIM, IK), \
66 : PROC//SK_": The condition `1 < size(sample, dim)` must hold. dim, shape(sample) = "//getStr([DIM, shape(sample, IK)]))
67 : #define CHECK_VAL_NDIM(PROC,DIM) \
68 : CHECK_ASSERTION(__LINE__, 0 < size(sample, DIM, IK), \
69 : PROC//SK_": The condition `0 < size(sample, dim)` must hold. dim, shape(sample) = "//getStr([DIM, shape(sample, IK)]))
70 : #define CHECK_VAL_DIM(PROC) \
71 : CHECK_ASSERTION(__LINE__, 1 <= dim .and. dim <= rank(sample), \
72 : PROC//SK_": The condition `1 <= dim .and. dim <= rank(sample)` must hold. dim, rank(sample) = "\
73 : //getStr([integer(IK) :: dim, rank(sample)]))
74 : #define CHECK_SUM_WEI(PROC) \
75 : CHECK_ASSERTION(__LINE__, 0 < sum(weight), \
76 : PROC//SK_": The condition `0 < sum(weight)` must hold. weight = "//getStr(weight))
77 : #define CHECK_VAL_WEI(PROC) \
78 : CHECK_ASSERTION(__LINE__, all(0._TKC <= weight), \
79 : PROC//SK_": The condition `all(0. <= weight)` must hold. weight = "//getStr(weight))
80 : #define CHECK_LEN_VAR(PROC) \
81 : CHECK_ASSERTION(__LINE__, size(sample, 3 - dim, IK) == size(var, 1, IK), \
82 : PROC//SK_": The condition `size(sample, 3 - dim) == size(var)` must hold. dim, size(sample, 3 - dim), size(var, 1) = "\
83 : //getStr([dim, size(sample, 3 - dim, IK), size(var, 1, IK)]))
84 : #define CHECK_VAL_MEANG(PROC,DIM) \
85 : CHECK_ASSERTION(__LINE__, all([minval(sample, DIM) <= meang .and. meang <= maxval(sample, DIM)]), \
86 : PROC//SK_": The condition `all([minval(sample, dim) <= meang .and. meang <= maxval(sample, dim)])` must hold. dim, minval(sample, dim), meang, maxval(sample, dim) = "//\
87 : getStr(DIM)//SK_"; "//getStr(minval(sample, DIM))//SK_"; "//getStr(meang)//SK_"; "//getStr(maxval(sample, DIM)))
88 : #define CHECK_LEN_MEANG(PROC) \
89 : CHECK_ASSERTION(__LINE__, size(sample, 3 - dim, IK) == size(meang, 1, IK), \
90 : PROC//SK_": The condition `size(sample, 3 - dim) == size(meang)` must hold. dim, size(sample, 3 - dim), size(meang, 1) = "\
91 : //getStr([dim, size(sample, 3 - dim, IK), size(meang, 1, IK)]))
92 : #define CHECK_LEN_MEAN(PROC) \
93 : CHECK_ASSERTION(__LINE__, size(sample, 3 - dim, IK) == size(mean, 1, IK), \
94 : PROC//SK_": The condition `size(sample, 3 - dim) == size(mean)` must hold. dim, size(sample, 3 - dim), size(mean, 1) = "\
95 : //getStr([dim, size(sample, 3 - dim, IK), size(mean, 1, IK)]))
96 : #define CHECK_LEN_WEI(PROC,DIM) \
97 : CHECK_ASSERTION(__LINE__, size(sample, DIM, IK) == size(weight, 1, IK), \
98 : PROC//SK_": The condition `size(sample, dim) == size(weight)` must hold. dim, size(sample, dim), size(weight, 1) = "\
99 : //getStr([DIM, size(sample, DIM, IK), size(weight, 1, IK)]))
100 : ! For very large sample, the following test is very crude and sensitive to the default vs. high-precision summation methods.
101 : #define CHECK_WEISUM(PROC) \
102 : CHECK_ASSERTION(__LINE__, abs(weisum - sum(weight)) < 1000 * epsilon(0._TKC), \
103 : PROC//SK_": The condition `0 < sum(weight)` must hold. weisum, sum(weight) = "//getStr([weisum, sum(weight)]))
104 :
105 : !%%%%%%%%%%%%%%%%%%%%%%%
106 : #if getVarCorrection_ENABLED
107 : !%%%%%%%%%%%%%%%%%%%%%%%
108 :
109 : #if Freq_ENABLED
110 9475 : CHECK_ASSERTION(__LINE__, 1._TKC < weisum, SK_"@setVarCorrection(): The condition `0 < weisum` must hold. weisum = "//getStr(weisum))
111 9475 : correction = weisum / (weisum - 1._TKC)
112 : #elif Reli_ENABLED
113 4855 : CHECK_ASSERTION(__LINE__, 0._TKC < weisum, SK_"@setVarCorrection(): The condition `0 < weisum` must hold. weisum = "//getStr(weisum))
114 4855 : CHECK_ASSERTION(__LINE__, 0._TKC < weisqs, SK_"@setVarCorrection(): The condition `0 < weisqs` must hold. weisqs = "//getStr(weisqs))
115 4855 : correction = weisum**2 / (weisum**2 - weisqs)
116 : #else
117 : #error "Unrecognized interface."
118 : #endif
119 :
120 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 : #elif getVarMerged_ENABLED && New_ENABLED
122 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 :
124 860 : call setVarMerged(varMerged, varB, varA, meanDiff, fracA)
125 :
126 : !%%%%%%%%%%%%%%%%%%%
127 : #elif setVarMerged_ENABLED
128 : !%%%%%%%%%%%%%%%%%%%
129 :
130 : ! Define the output value for setVarMerged.
131 : #if Old_ENABLED
132 : #define CHECK_LEN_TARGET_VAR
133 : #define TARGET_VAR varB
134 : #elif New_ENABLED
135 : #define TARGET_VAR varMerged
136 : #define CHECK_LEN_TARGET_VAR \
137 : CHECK_ASSERTION(__LINE__, size(varMerged, 1, IK) == size(varB, 1, IK), SK_"@setVarMerged(): The condition `size(varMerged) == size(varB)` must hold. size(varMerged), size(varB) = "//getStr([size(varMerged, 1, IK), size(varB, 1, IK)]))
138 : #else
139 : #error "Unrecognized interface."
140 : #endif
141 : #if D1_ENABLED
142 : integer(IK) :: idim
143 : real(TKC) :: fracAB
144 : #endif
145 : real(TKC) :: fracB
146 2580 : fracB = 1._TKC - fracA
147 2580 : CHECK_ASSERTION(__LINE__, 0._TKC < fracA .and. fracA < 1._TKC, SK_"@setVarMerged(): The condition `0 < fracA .and. fracA < 1` must hold. fracA = "//getStr(fracA))
148 : #if D0_ENABLED
149 1290 : TARGET_VAR = fracA * varA + fracB * varB + fracA * fracB * GET_ABSQ(meanDiff)
150 : #elif D1_ENABLED
151 2580 : CHECK_LEN_TARGET_VAR
152 3870 : CHECK_ASSERTION(__LINE__, size(varB, 1, IK) == size(varA, 1, IK), SK_"@setVarMerged(): The condition `size(varB) == size(varA)` must hold. size(varB), size(varA) = "//getStr([size(varB, 1, IK), size(varA, 1, IK)]))
153 3870 : CHECK_ASSERTION(__LINE__, size(meanDiff, 1, IK) == size(varA, 1, IK), SK_"@setVarMerged(): The condition `size(meanDiff) == size(varA)` must hold. size(meanDiff), size(varA) = "//getStr([size(meanDiff, 1, IK), size(varA, 1, IK)]))
154 1290 : fracAB = fracA * fracB
155 : do concurrent(idim = 1 : size(varA, 1, IK))
156 5203 : TARGET_VAR(idim) = fracB * varB(idim) + fracA * varA(idim) + fracAB * GET_ABSQ(meanDiff(idim))
157 : end do
158 : #else
159 : #error "Unrecognized interface."
160 : #endif
161 : #undef CHECK_LEN_TARGET_VAR
162 : #undef TARGET_VAR
163 :
164 : !%%%%%%%%%%%%%%%%%%%%%%%
165 : #elif setVarMeanMerged_ENABLED
166 : !%%%%%%%%%%%%%%%%%%%%%%%
167 :
168 : ! Define the output value for setVarMerged.
169 : #if Old_ENABLED
170 : #define CHECK_LEN_TARGET_AVG
171 : #define CHECK_LEN_TARGET_VAR
172 : #define TARGET_AVG meanB
173 : #define TARGET_VAR varB
174 : #elif New_ENABLED
175 : #define CHECK_LEN_TARGET_AVG \
176 : CHECK_ASSERTION(__LINE__, size(meanMerged, 1, IK) == size(varB, 1, IK), SK_"@setVarMeanMerged(): The condition `size(meanMerged) == size(varB)` must hold. size(meanMerged), size(varB) = "//getStr([size(meanMerged, 1, IK), size(varB, 1, IK)]))
177 : #define CHECK_LEN_TARGET_VAR \
178 : CHECK_ASSERTION(__LINE__, size(varMerged, 1, IK) == size(varB, 1, IK), SK_"@setVarMeanMerged(): The condition `size(varMerged) == size(varB)` must hold. size(varMerged), size(varB) = "//getStr([size(varMerged, 1, IK), size(varB, 1, IK)]))
179 : #define TARGET_AVG meanMerged
180 : #define TARGET_VAR varMerged
181 : #else
182 : #error "Unrecognized interface."
183 : #endif
184 : TYPE_OF_SAMPLE :: temp
185 : #if D1_ENABLED
186 : integer(IK) :: idim
187 : real(TKC) :: fracAB
188 : #endif
189 : real(TKC) :: fracB
190 1720 : fracB = 1._TKC - fracA
191 1720 : CHECK_ASSERTION(__LINE__, 0._TKC < fracA .and. fracA < 1._TKC, SK_"@setVarMeanMerged(): The condition `0 < fracA .and. fracA < 1` must hold. fracA = "//getStr(fracA))
192 : #if D0_ENABLED
193 860 : temp = meanA - meanB
194 860 : TARGET_VAR = fracA * varA + fracB * varB + fracA * fracB * GET_ABSQ(temp)
195 : #elif D1_ENABLED
196 1290 : CHECK_LEN_TARGET_VAR
197 2580 : CHECK_ASSERTION(__LINE__, size(varB, 1, IK) == size(varA, 1, IK), SK_"@setVarMeanMerged(): The condition `size(varB) == size(varA)` must hold. size(varB), size(varA) = "//getStr([size(varB, 1, IK), size(varA, 1, IK)]))
198 2580 : CHECK_ASSERTION(__LINE__, size(meanB, 1, IK) == size(varA, 1, IK), SK_"@setVarMeanMerged(): The condition `size(meanB) == size(varA)` must hold. size(meanB), size(varA) = "//getStr([size(meanB, 1, IK), size(varA, 1, IK)]))
199 2580 : CHECK_ASSERTION(__LINE__, size(meanA, 1, IK) == size(varA, 1, IK), SK_"@setVarMeanMerged(): The condition `size(meanA) == size(varA)` must hold. size(meanA), size(varA) = "//getStr([size(meanA, 1, IK), size(varA, 1, IK)]))
200 860 : fracAB = fracA * fracB
201 : do concurrent(idim = 1 : size(varA, 1, IK))
202 2558 : temp = meanA(idim) - meanB(idim)
203 3418 : TARGET_VAR(idim) = fracB * varB(idim) + fracA * varA(idim) + fracAB * GET_ABSQ(temp)
204 : end do
205 : #else
206 : #error "Unrecognized interface."
207 : #endif
208 3418 : TARGET_AVG = fracA * meanA + fracB * meanB
209 : #undef CHECK_LEN_TARGET_AVG
210 : #undef CHECK_LEN_TARGET_VAR
211 : #undef TARGET_AVG
212 : #undef TARGET_VAR
213 :
214 : !%%%%%%%%%%%%%
215 : #elif getVar_ENABLED
216 : !%%%%%%%%%%%%%
217 :
218 : ! Define the default dimension.
219 : #if ALL_ENABLED
220 : #define DIM_ARG
221 : #elif DIM_ENABLED
222 : #define DIM_ARG, dim
223 : #else
224 : #error "Unrecognized interface."
225 : #endif
226 : real(TKC) :: normfac
227 : ! Define the weight sum.
228 : #if WNO_ENABLED
229 : real(TKC) :: weisum
230 34432 : weisum = real(product(shape(sample, IK)), TKC)
231 12337 : call setVar(var, getMean(sample DIM_ARG), sample DIM_ARG)
232 : #elif WTI_ENABLED || WTR_ENABLED
233 : TYPE_OF_WEIGHT :: weisum
234 75663 : weisum = sum(weight)
235 10003 : call setVar(var, getMean(sample DIM_ARG, weight), sample DIM_ARG, weight, weisum)
236 : #else
237 : #error "Unrecognized interface."
238 : #endif
239 22340 : if (present(correction)) then
240 3630 : CHECK_ASSERTION(__LINE__, same_type_as(correction, fweight) .or. same_type_as(correction, rweight), SK_"@getVar(): The condition `same_type_as(correction, fweight) .or. same_type_as(correction, rweight)` must hold.")
241 : #if WNO_ENABLED
242 420 : normfac = getVarCorrection(real(weisum, TKC))
243 : #elif (WTI_ENABLED || WTR_ENABLED)
244 3210 : if (same_type_as(correction, fweight)) then
245 1605 : normfac = getVarCorrection(real(weisum, TKC))
246 1605 : elseif (same_type_as(correction, rweight)) then
247 12128 : normfac = getVarCorrection(real(weisum, TKC), real(sum(weight**2), TKC))
248 : end if
249 : #else
250 : #error "Unrecognized interface."
251 : #endif
252 7257 : var = var * normfac
253 : end if
254 : #undef DIM_ARG
255 :
256 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
257 : #elif setVar_ENABLED && D1_ENABLED
258 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
259 :
260 : TYPE_OF_SAMPLE :: temp
261 : integer(IK) :: isam
262 : #if WNO_ENABLED
263 : real(TKC) :: weisum
264 4987 : weisum = real(size(sample, 1, IK), TKC)
265 : #elif WTI_ENABLED || WTR_ENABLED
266 117242 : CHECK_WEISUM(SK_"@setVar()")
267 53722 : CHECK_SUM_WEI(SK_"@setVar()")
268 53722 : CHECK_VAL_WEI(SK_"@setVar()")
269 39192 : CHECK_LEN_WEI(SK_"@setVar()",1_IK)
270 : #else
271 : #error "Unrecognized interface."
272 : #endif
273 59140 : CHECK_VAL_NSAM(SK_"@setVar()",1_IK)
274 14785 : var = 0._TKC
275 2261736 : do isam = 1, size(sample, 1, IK)
276 2246951 : temp = GET_SHIFTED(sample(isam),mean)
277 2261736 : var = var + GET_WEIGHTED(GET_ABSQ(temp),weight(isam))
278 : end do
279 14785 : var = var / weisum
280 :
281 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282 : #elif setVar_ENABLED && D2_ENABLED && ALL_ENABLED
283 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
284 :
285 : TYPE_OF_SAMPLE :: temp
286 : integer(IK) :: idim, jdim
287 : integer(IK) :: mdim, ndim
288 : #if WNO_ENABLED
289 : real(TKC) :: weisum
290 2412 : weisum = size(sample, kind = IK)
291 : #elif WTI_ENABLED || WTR_ENABLED
292 : integer(IK) :: iwei
293 : iwei = 0_IK
294 28830 : CHECK_ASSERTION(__LINE__, size(sample, kind = IK) == size(weight, 1, IK), SK_"@setVar(): The condition `size(sample) == size(weight)` must hold. shape(sample), size(weight) = "//getStr([shape(sample, IK), size(weight, 1, IK)]))
295 68771 : CHECK_SUM_WEI(SK_"@setVar()")
296 68771 : CHECK_VAL_WEI(SK_"@setVar()")
297 142347 : CHECK_WEISUM(SK_"@setVar()")
298 : #else
299 : #error "Unrecognized interface."
300 : #endif
301 7217 : CHECK_ASSERTION(__LINE__, 1_IK < size(sample, kind = IK), SK_"@setVar(): The condition `1 < size(sample)` must hold. shape(sample) = "//getStr(shape(sample, IK)))
302 7217 : mdim = size(sample, 1, IK)
303 7217 : ndim = size(sample, 2, IK)
304 7217 : var = 0._TKC
305 34329 : do jdim = 1, ndim
306 130468 : do idim = 1, mdim
307 : #if WTI_ENABLED || WTR_ENABLED
308 63966 : iwei = iwei + 1_IK
309 : #endif
310 96139 : temp = GET_SHIFTED(sample(idim, jdim),mean)
311 123251 : var = var + GET_WEIGHTED(GET_ABSQ(temp),weight(iwei))
312 : end do
313 : end do
314 7217 : var = var / weisum
315 :
316 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
317 : #elif setVar_ENABLED && D2_ENABLED && DIM_ENABLED
318 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
319 :
320 : real(TKC) :: normFac
321 : TYPE_OF_SAMPLE :: temp
322 : integer(IK) :: idim, isam, ndim, nsam
323 14766 : ndim = size(sample, 3 - dim, IK)
324 5008 : nsam = size(sample, dim, IK)
325 : #if WNO_ENABLED
326 9758 : normFac = 1._TKC / real(nsam, TKC)
327 : #elif WTI_ENABLED || WTR_ENABLED
328 60348 : CHECK_WEISUM(SK_"@setVar()")
329 27670 : CHECK_SUM_WEI(SK_"@setVar()")
330 27670 : CHECK_VAL_WEI(SK_"@setVar()")
331 20032 : CHECK_LEN_WEI(SK_"@setVar()",dim)
332 5008 : normFac = 1._TKC / real(weisum, TKC)
333 : #else
334 : #error "Unrecognized interface."
335 : #endif
336 44298 : CHECK_VAL_DIM(SK_"@setVar()")
337 59064 : CHECK_LEN_VAR(SK_"@setVar()")
338 88596 : CHECK_VAL_NSAM(SK_"@setVar()",dim)
339 88596 : CHECK_VAL_NDIM(SK_"@setVar()",3 - dim)
340 : #if Avg_ENABLED
341 54248 : CHECK_LEN_MEAN(SK_"@setVar()")
342 : #endif
343 14766 : if (dim == 2_IK) then
344 37810 : var = 0._TKC
345 64124 : do isam = 1, nsam
346 11291 : do concurrent(idim = 1 : ndim)
347 124048 : temp = GET_SHIFTED(sample(idim, isam),mean(idim))
348 176881 : var(idim) = var(idim) + GET_WEIGHTED(GET_ABSQ(temp),weight(isam))
349 : end do
350 : end do
351 37810 : var = var * normFac
352 : else ! dim = 1
353 13890 : do idim = 1, ndim
354 10415 : var(idim) = 0._TKC
355 56862 : do isam = 1, nsam
356 46447 : temp = GET_SHIFTED(sample(isam, idim),mean(idim))
357 56862 : var(idim) = var(idim) + GET_WEIGHTED(GET_ABSQ(temp),weight(isam))
358 : end do
359 13890 : var(idim) = var(idim) * normFac
360 : end do
361 : end if ! dim
362 :
363 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
364 : #elif setVarMean_ENABLED && D1_ENABLED
365 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
366 :
367 : integer(IK) :: isam, nsam
368 : TYPE_OF_SAMPLE :: temp1
369 : #if WNO_ENABLED
370 : real(TKC) :: weisum
371 402 : weisum = size(sample, 1, IK)
372 : #elif WTI_ENABLED || WTR_ENABLED
373 : TYPE_OF_SAMPLE :: temp2
374 6416 : CHECK_LEN_WEI(SK_"@setVarMean()",1_IK)
375 8704 : CHECK_SUM_WEI(SK_"@setVarMean()")
376 8704 : CHECK_VAL_WEI(SK_"@setVarMean()")
377 802 : weisum = 0
378 : #else
379 : #error "Unrecognized interface."
380 : #endif
381 29724 : CHECK_VAL_MEANG(SK_"@setVarMean()",1_IK)
382 9632 : CHECK_VAL_NSAM(SK_"@setVarMean()",1_IK)
383 : nsam = size(sample, 1, IK)
384 2408 : mean = 0._TKC
385 2408 : var = 0._TKC
386 13068 : do isam = 1, nsam
387 10660 : temp1 = sample(isam) - meang
388 : #if WNO_ENABLED
389 3560 : var = var + GET_ABSQ(temp1)
390 4364 : mean = mean + temp1
391 : #elif WTI_ENABLED || WTR_ENABLED
392 7100 : temp2 = temp1 * weight(isam)
393 7100 : weisum = weisum + weight(isam)
394 7100 : var = var + GET_PROD(temp1,temp2)
395 8704 : mean = mean + temp2
396 : #endif
397 : end do
398 2408 : mean = mean / weisum
399 2408 : var = (var - GET_ABSQ(mean) * weisum) / weisum
400 2408 : mean = mean + meang
401 :
402 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
403 : #elif setVarMean_ENABLED && D2_ENABLED && ALL_ENABLED
404 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
405 :
406 : integer(IK) :: idim, jdim
407 : integer(IK) :: mdim, ndim
408 : TYPE_OF_SAMPLE :: temp1
409 : #if WNO_ENABLED
410 : real(TKC) :: weisum
411 402 : weisum = size(sample, kind = IK)
412 : #elif WTI_ENABLED || WTR_ENABLED
413 : integer(IK) :: iwei
414 : TYPE_OF_SAMPLE :: temp2
415 4812 : CHECK_ASSERTION(__LINE__, size(sample, kind = IK) == size(weight, 1, IK), SK_"@setVarMean(): The condition `size(sample) == size(weight)` must hold. shape(sample), size(weight) = "//getStr([shape(sample, IK), size(weight, 1, IK)]))
416 11434 : CHECK_SUM_WEI(SK_"@setVarMean()")
417 11434 : CHECK_VAL_WEI(SK_"@setVarMean()")
418 : iwei = 0_IK
419 802 : weisum = 0
420 : #else
421 : #error "Unrecognized interface."
422 : #endif
423 1204 : CHECK_ASSERTION(__LINE__, 1_IK < size(sample, kind = IK), SK_"@setVarMean(): The condition `1 < size(sample)` must hold. shape(sample) = "//getStr(shape(sample, IK)))
424 46195 : CHECK_ASSERTION(__LINE__, minval(sample) <= meang .and. meang <= maxval(sample), SK_"@setVarMean(): The condition `minval(sample) <= meang .and. meang <= minval(sample)` must hold. minval(sample), meang, maxval(sample) = "//getStr([minval(sample), meang, maxval(sample)]))
425 1204 : mdim = size(sample, 1, IK)
426 1204 : ndim = size(sample, 2, IK)
427 1204 : mean = 0._TKC
428 1204 : var = 0._TKC
429 5601 : do jdim = 1, ndim
430 21569 : do idim = 1, mdim
431 15968 : temp1 = sample(idim, jdim) - meang
432 : #if WNO_ENABLED
433 5336 : var = var + GET_ABSQ(temp1)
434 6805 : mean = mean + temp1
435 : #elif WTI_ENABLED || WTR_ENABLED
436 10632 : iwei = iwei + 1_IK
437 10632 : temp2 = temp1 * weight(iwei)
438 10632 : weisum = weisum + weight(iwei)
439 10632 : var = var + GET_PROD(temp1,temp2)
440 13560 : mean = mean + temp2
441 : #endif
442 : end do
443 : end do
444 1204 : mean = mean / weisum
445 1204 : var = (var - GET_ABSQ(mean) * weisum) / weisum
446 1204 : mean = mean + meang
447 :
448 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
449 : #elif setVarMean_ENABLED && D2_ENABLED && DIM_ENABLED
450 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
451 :
452 : real(TKC) :: normFac
453 : TYPE_OF_SAMPLE :: temp1
454 : #if WTI_ENABLED || WTR_ENABLED
455 : TYPE_OF_SAMPLE :: temp2
456 : #endif
457 : integer(IK) :: idim, isam, ndim, nsam
458 1206 : ndim = size(sample, 3 - dim, IK)
459 1206 : nsam = size(sample, dim, IK)
460 3618 : CHECK_VAL_DIM(SK_"@setVarMean()")
461 4824 : CHECK_LEN_VAR(SK_"@setVarMean()")
462 7236 : CHECK_VAL_NSAM(SK_"@setVarMean()",dim)
463 7236 : CHECK_VAL_NDIM(SK_"@setVarMean()",3 - dim)
464 4884 : CHECK_VAL_MEANG(SK_"@setVarMean()",dim)
465 4824 : CHECK_LEN_MEAN(SK_"@setVarMean()")
466 : #if WNO_ENABLED
467 404 : normFac = 1._TKC / real(nsam, TKC)
468 404 : if (dim == 2_IK) then
469 809 : var = 0._TKC
470 809 : mean = 0._TKC
471 1062 : do isam = 1, nsam
472 3777 : do idim = 1, ndim
473 2715 : temp1 = sample(idim, isam) - meang(idim)
474 2715 : var(idim) = var(idim) + GET_ABSQ(temp1)
475 3581 : mean(idim) = mean(idim) + temp1
476 : end do
477 : end do
478 : do concurrent(idim = 1 : ndim)
479 613 : mean(idim) = mean(idim) * normFac
480 613 : var(idim) = (var(idim) - GET_ABSQ(mean(idim)) * nsam) * normFac
481 809 : mean(idim) = mean(idim) + meang(idim)
482 : end do
483 : else
484 : do concurrent(idim = 1 : ndim)
485 622 : var(idim) = 0._TKC
486 622 : mean(idim) = 0._TKC
487 3393 : do isam = 1, nsam
488 2771 : temp1 = sample(isam, idim) - meang(idim)
489 2771 : var(idim) = var(idim) + GET_ABSQ(temp1)
490 3393 : mean(idim) = mean(idim) + temp1
491 : end do
492 622 : mean(idim) = mean(idim) * normFac
493 622 : var(idim) = (var(idim) - GET_ABSQ(mean(idim)) * nsam) * normFac
494 830 : mean(idim) = mean(idim) + meang(idim)
495 : end do
496 : end if
497 : #elif WTI_ENABLED || WTR_ENABLED
498 3208 : CHECK_LEN_WEI(SK_"@setVarMean()",dim)
499 4351 : CHECK_SUM_WEI(SK_"@setVarMean()")
500 4351 : CHECK_VAL_WEI(SK_"@setVarMean()")
501 802 : weisum = 0
502 802 : if (dim == 2_IK) then
503 1603 : var = 0._TKC
504 1603 : mean = 0._TKC
505 2106 : do isam = 1, nsam
506 1717 : weisum = weisum + weight(isam)
507 389 : do concurrent(idim = 1 : ndim)
508 5370 : temp1 = sample(idim, isam) - meang(idim)
509 5370 : temp2 = temp1 * weight(isam)
510 5370 : var(idim) = var(idim) + GET_PROD(temp1,temp2)
511 7087 : mean(idim) = mean(idim) + temp2
512 : end do
513 : end do
514 389 : normFac = 1._TKC / weisum
515 : do concurrent(idim = 1 : ndim)
516 1214 : mean(idim) = mean(idim) * normFac
517 1214 : var(idim) = (var(idim) - GET_ABSQ(mean(idim)) * weisum) * normFac
518 1603 : mean(idim) = mean(idim) + meang(idim)
519 : end do
520 : else
521 : idim = 1_IK
522 413 : var(idim) = 0._TKC
523 413 : mean(idim) = 0._TKC
524 2245 : do isam = 1, nsam
525 1832 : weisum = weisum + weight(isam)
526 1832 : temp1 = sample(isam, idim) - meang(idim)
527 1832 : temp2 = temp1 * weight(isam)
528 1832 : var(idim) = var(idim) + GET_PROD(temp1,temp2)
529 2245 : mean(idim) = mean(idim) + temp2
530 : end do
531 413 : normFac = 1._TKC / real(weisum, TKC)
532 413 : mean(idim) = mean(idim) * normFac
533 413 : var(idim) = (var(idim) - GET_ABSQ(mean(idim)) * weisum) * normFac
534 413 : mean(idim) = mean(idim) + meang(idim)
535 1229 : do idim = 2, ndim
536 816 : var(idim) = 0._TKC
537 816 : mean(idim) = 0._TKC
538 4466 : do isam = 1, nsam
539 3650 : temp1 = sample(isam, idim) - meang(idim)
540 3650 : temp2 = temp1 * weight(isam)
541 3650 : var(idim) = var(idim) + GET_PROD(temp1,temp2)
542 4466 : mean(idim) = mean(idim) + temp2
543 : end do
544 816 : mean(idim) = mean(idim) * normFac
545 816 : var(idim) = (var(idim) - GET_ABSQ(mean(idim)) * weisum) * normFac
546 1229 : mean(idim) = mean(idim) + meang(idim)
547 : end do
548 : end if
549 : #else
550 : #error "Unrecognized interface."
551 : #endif
552 :
553 : #else
554 : !%%%%%%%%%%%%%%%%%%%%%%%%
555 : #error "Unrecognized interface."
556 : !%%%%%%%%%%%%%%%%%%%%%%%%
557 : #endif
558 : #undef CHECK_VAL_MEANG
559 : #undef CHECK_LEN_MEANG
560 : #undef CHECK_LEN_MEAN
561 : #undef TYPE_OF_SAMPLE
562 : #undef TYPE_OF_WEIGHT
563 : #undef CHECK_VAL_NSAM
564 : #undef CHECK_VAL_NDIM
565 : #undef CHECK_LEN_WEI
566 : #undef CHECK_VAL_DIM
567 : #undef CHECK_LEN_VAR
568 : #undef CHECK_SUM_WEI
569 : #undef CHECK_VAL_WEI
570 : #undef CHECK_WEISUM
571 : #undef GET_WEIGHTED
572 : #undef GET_SHIFTED
573 : #undef GET_ABSQ
574 : #undef GET_PROD
575 : #undef FIRST
|
