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_sampleCCF](@ref pm_sampleCCF).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Saturday 4:40 PM, August 21, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the type and kind and the conjugation rule.
28 : #if CK_ENABLED
29 : #define GET_ABSQ(X)(real(X)**2 + aimag(X)**2)
30 : #define TYPE_OF_SEQ complex(TKC)
31 : #define GET_RE(X)X%re
32 : #elif RK_ENABLED
33 : #define TYPE_OF_SEQ real(TKC)
34 : #define GET_ABSQ(X)X**2
35 : #define GET_RE(X)X
36 : #else
37 : #error "Unrecognized interface."
38 : #endif
39 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 : #if getACF_ENABLED && D1_ENABLED
41 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 :
43 : real(TKC) :: normfac
44 : TYPE_OF_SEQ :: meanf
45 16016 : class(*), allocatable :: norm_def
46 : TYPE_OF_SEQ, parameter :: ZERO = 0._TKC
47 : TYPE_OF_SEQ, allocatable :: ff(:), coef(:)
48 : integer(IK), allocatable :: factor(:)
49 : logical(LK) :: unscaled, inf
50 : integer(IK) :: lagmax
51 : integer(IK) :: lagmin
52 : integer(IK) :: lenacf
53 :
54 9608 : lagmax = size(f, 1, IK) - 1_IK
55 9608 : lagmin = 1_IK - size(f, 1, IK)
56 9608 : if (present(lag)) then
57 6402 : if (size(lag, 1, IK) == 0_IK) then
58 3200 : allocate(acf(0))
59 : return
60 : end if
61 3202 : CHECK_ASSERTION(__LINE__, isAscending(lag), SK_"@setACF(): The condition `isAscending(lag)` must hold. lag = "//getStr(lag))
62 3202 : lagmax = max(lagmax, lag(size(lag, 1, IK)))
63 3202 : lagmin = min(lagmin, lag(1))
64 : end if
65 6408 : lenacf = lagmax - lagmin + 1_IK
66 :
67 : ! Allocate and pad the arrays.
68 :
69 6408 : allocate(ff(lenacf), coef(lenacf), acf(lenacf))
70 :
71 : ! Shift the padded arrays if needed.
72 :
73 6408 : if (present(norm)) then
74 4804 : unscaled = logical(same_type_as(norm, zscore_type()) .or. same_type_as(norm, stdscale_type()), LK)
75 4804 : norm_def = norm
76 : else
77 : unscaled = .true._LK
78 1604 : norm_def = zscore_type() ! Intel cannot digest `same_type_as(norm_def, zscore_type())` if norm_def is set to constant.
79 : end if
80 :
81 6408 : if (same_type_as(norm_def, meanshift_type()) .or. same_type_as(norm_def, zscore_type())) then
82 4804 : meanf = getMean(f)
83 105966 : ff(1 : size(f, 1, IK)) = f - meanf
84 1604 : elseif (same_type_as(norm_def, stdscale_type()) .or. same_type_as(norm_def, nothing_type())) then
85 72258 : ff(1 : size(f, 1, IK)) = f
86 : else
87 : error stop MODULE_NAME//SK_"@getACF(): Unrecognized `norm`." ! LCOV_EXCL_LINE
88 : end if
89 267684 : ff(size(f, 1, IK) + 1 : lenacf) = ZERO
90 :
91 : ! Compute the ACF.
92 :
93 25080 : factor = getFactorFFT(ff, coef)
94 6408 : call setACF(factor, coef, ff, acf, inf)
95 :
96 : ! Normalize and shift.
97 :
98 6408 : normfac = 1._TKC / lenacf ! default scale.
99 6408 : if (inf) then
100 2468 : if (unscaled) normfac = 1._TKC / GET_RE(ff(1))
101 2468 : if (present(lag)) then
102 106108 : acf = cshift(ff * normfac, lagmin)
103 : else
104 23210 : acf = ff(1 : 1 + lagmax) * normfac
105 : end if
106 : else
107 3940 : if (unscaled) normfac = 1._TKC / GET_RE(acf(1))
108 3940 : if (present(lag)) then
109 124588 : acf = cshift(acf * normfac, lagmin)
110 : else
111 170776 : acf = acf(1 : 1 + lagmax) * normfac
112 : end if
113 : end if
114 490096 : if (present(lag)) acf = acf(lag - lagmin + 1)
115 6408 : deallocate(ff, coef) ! essential for gfortran on heap.
116 :
117 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118 : #elif setACF_ENABLED && FP_ENABLED && D1_ENABLED
119 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120 :
121 : integer(IK) :: isam, nsam
122 6577 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) >= 2_IK, SK_"@setACF(): The condition `1 < size(f)` must hold. size(f) = "//getStr(size(f, 1, IK)))
123 19731 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) == size(coef, 1, IK), SK_"@setACF(): The condition `size(f) == size(coef)` must hold. size(f), size(coef) = "//getStr([size(f, 1, IK), size(coef, 1, IK)]))
124 19731 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) == size(work, 1, IK), SK_"@setACF(): The condition `size(f) == size(work)` must hold. size(f), size(work) = "//getStr([size(f, 1, IK), size(work, 1, IK)]))
125 6577 : call setFFTF(factor, coef, f, work, inf)
126 : nsam = size(f, 1, IK)
127 : #if CK_ENABLED
128 3286 : if (inf) then
129 : ! fft of `f` is in `work`.
130 : do concurrent(isam = 1 : nsam)
131 108492 : f(isam)%re = work(isam)%re**2 + work(isam)%im**2
132 110536 : f(isam)%im = 0._TKC
133 : end do
134 : else
135 : ! fft of `f` is in `f`.
136 : do concurrent(isam = 1 : nsam)
137 73822 : f(isam)%re = f(isam)%re**2 + f(isam)%im**2
138 75064 : f(isam)%im = 0._TKC
139 : end do
140 : end if
141 : #elif RK_ENABLED
142 3291 : if (mod(nsam, 2_IK) == 0_IK) nsam = nsam - 1_IK
143 3291 : if (inf) then
144 : ! fft of `f` is in `work`.
145 2024 : f(1) = work(1)**2
146 : do concurrent(isam = 2 : nsam : 2)
147 90876 : f(isam) = work(isam)**2 + work(isam + 1)**2
148 92900 : f(isam + 1) = 0._TKC
149 : end do
150 2024 : if (nsam < size(f, 1, IK)) f(nsam + 1) = work(nsam + 1)**2
151 : else
152 : ! fft of `f` is in `f`.
153 1267 : f(1) = f(1)**2
154 : do concurrent(isam = 2 : nsam : 2)
155 36092 : f(isam) = f(isam)**2 + f(isam + 1)**2
156 37359 : f(isam + 1) = 0._TKC
157 : end do
158 1267 : if (nsam < size(f, 1, IK)) f(nsam + 1) = f(nsam + 1)**2
159 : end if
160 : #else
161 : #error "Unrecognized interface."
162 : #endif
163 : ! result is now in `f`.
164 6577 : call setFFTR(factor, coef, f, work, inf)
165 6577 : inf = .not. inf
166 : #undef SET_CCF
167 :
168 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
169 : #elif getCCF_ENABLED && FG_ENABLED
170 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%
171 :
172 : TYPE_OF_SEQ :: meanf, meang
173 : TYPE_OF_SEQ, parameter :: ZERO = 0._TKC
174 : TYPE_OF_SEQ, allocatable :: ff(:), gg(:), coef(:)
175 : real(TKC) :: normfac, sumsqf, sumsqg
176 : integer(IK), allocatable :: factor(:)
177 32812 : class(*), allocatable :: norm_def
178 : integer(IK) :: lagmax
179 : integer(IK) :: lagmin
180 : integer(IK) :: lenccf
181 : integer(IK) :: iell
182 : logical(LK) :: inf
183 :
184 18406 : lagmax = size(g, 1, IK) - 1_IK
185 18406 : lagmin = 1_IK - size(f, 1, IK)
186 18406 : if (present(lag)) then
187 14404 : if (size(lag, 1, IK) == 0_IK) then
188 4000 : allocate(ccf(0))
189 : return
190 : end if
191 10404 : CHECK_ASSERTION(__LINE__, isAscending(lag), SK_"@setCCF(): The condition `isAscending(lag)` must hold. lag = "//getStr(lag))
192 10404 : lagmax = max(lagmax, lag(size(lag, 1, IK)))
193 10404 : lagmin = min(lagmin, lag(1))
194 : end if
195 14406 : lenccf = lagmax - lagmin + 1_IK
196 :
197 : ! Allocate and pad the arrays.
198 :
199 14406 : allocate(ff(lenccf), gg(lenccf), coef(lenccf), ccf(lenccf))
200 :
201 : ! Shift the padded arrays if needed.
202 :
203 14406 : if (present(norm)) then
204 11201 : norm_def = norm
205 : else
206 3205 : norm_def = zscore_type()
207 : end if
208 :
209 14406 : normfac = 1._TKC / lenccf
210 14406 : if (same_type_as(norm_def, meanshift_type()) .or. same_type_as(norm_def, zscore_type())) then
211 9605 : meanf = getMean(f)
212 9605 : meang = getMean(g)
213 212724 : ff(1 : size(f, 1, IK)) = f - meanf
214 211086 : gg(1 : size(g, 1, IK)) = g - meang
215 : else
216 106452 : ff(1 : size(f, 1, IK)) = f
217 105360 : gg(1 : size(g, 1, IK)) = g
218 : end if
219 518868 : ff(size(f, 1, IK) + 1 : lenccf) = ZERO
220 521598 : gg(size(g, 1, IK) + 1 : lenccf) = ZERO
221 :
222 14406 : if (same_type_as(norm_def, stdscale_type()) .or. same_type_as(norm_def, zscore_type())) then
223 4800 : sumsqf = 0._TKC
224 212766 : do iell = 1, size(f, 1, IK)
225 212766 : sumsqf = sumsqf + GET_ABSQ(ff(iell))
226 : end do
227 4800 : sumsqg = 0._TKC
228 211128 : do iell = 1, size(g, 1, IK)
229 211128 : sumsqg = sumsqg + GET_ABSQ(gg(iell))
230 : end do
231 9606 : normfac = normfac / sqrt(sumsqf * sumsqg)
232 4800 : elseif (.not. (same_type_as(norm_def, meanshift_type()) .or. same_type_as(norm_def, nothing_type()))) then
233 : error stop MODULE_NAME//SK_"@getCCF(): Unrecognized `norm`." ! LCOV_EXCL_LINE
234 : end if
235 :
236 58630 : factor = getFactorFFT(ff, coef)
237 14406 : call setCCF(factor, coef, ff, gg, ccf, inf)
238 :
239 14406 : if (inf) then
240 360712 : ccf = cshift(ff * normfac, lagmin)
241 : else
242 462926 : ccf = cshift(gg * normfac, lagmin)
243 : end if
244 1307538 : if (present(lag)) ccf = ccf(lag - lagmin + 1)
245 14406 : deallocate(ff, gg, coef) ! essential for gfortran on heap.
246 :
247 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
248 : #elif setCCF_ENABLED && FP_ENABLED && FG_ENABLED
249 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
250 :
251 : integer(IK) :: isam, nsam
252 22744 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) >= 2_IK, SK_"@setCCF(): The condition `1 < size(f)` must hold. size(f) = "//getStr(size(f, 1, IK)))
253 68232 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) == size(g, 1, IK), SK_"@setCCF(): The condition `size(f) == size(g)` must hold. size(f), size(g) = "//getStr([size(f, 1, IK), size(g, 1, IK)]))
254 68232 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) == size(coef, 1, IK), SK_"@setCCF(): The condition `size(f) == size(coef)` must hold. size(f), size(coef) = "//getStr([size(f, 1, IK), size(coef, 1, IK)]))
255 68232 : CHECK_ASSERTION(__LINE__, size(f, 1, IK) == size(work, 1, IK), SK_"@setCCF(): The condition `size(f) == size(work)` must hold. size(f), size(work) = "//getStr([size(f, 1, IK), size(work, 1, IK)]))
256 22744 : call setFFTF(factor, coef, f, work, inf)
257 : nsam = size(f, 1, IK)
258 : #if CK_ENABLED
259 : #define SET_CCF(RES,TS1,TS2)\
260 : do concurrent(isam = 1 : nsam); \
261 : RES(isam) = TS1(isam) * conjg(TS2(isam)); \
262 : end do;
263 11370 : if (inf) then
264 : ! fft of `f` is in `work`.
265 6458 : call setFFTF(factor, coef, g, f, inf)
266 6458 : if (inf) then
267 : ! fft of `g` is in `f`.
268 357102 : SET_CCF(f,work,f)
269 : else
270 : ! fft of `g` is in `g`.
271 0 : SET_CCF(f,work,g)
272 : end if
273 : else
274 : ! fft of `f` is in `f`.
275 4912 : call setFFTF(factor, coef, g, work, inf)
276 4912 : if (inf) then
277 : ! fft of `g` is in `work`.
278 0 : SET_CCF(f,f,work)
279 : else
280 : ! fft of `g` is in `g`.
281 287298 : SET_CCF(f,f,g)
282 : end if
283 : end if
284 : #elif RK_ENABLED
285 : #define SET_CCF(RES,TS1,TS2)\
286 : RES(1) = TS2(1) * TS1(1); \
287 : do concurrent(isam = 2 : nsam : 2); \
288 : t1 = TS1(isam) * TS2(isam) + TS1(isam + 1) * TS2(isam + 1); \
289 : t2 = TS1(isam + 1) * TS2(isam) - TS1(isam) * TS2(isam + 1); \
290 : RES(isam + 1) = t2; \
291 : RES(isam) = t1; \
292 : end do; \
293 : if (nsam < size(f, 1, IK)) f(nsam + 1) = TS2(nsam + 1) * TS1(nsam + 1);
294 : block
295 : TYPE_OF_SEQ :: t1, t2
296 11374 : if (mod(nsam, 2_IK) == 0_IK) nsam = nsam - 1_IK
297 11374 : if (inf) then
298 : ! fft of `f` is in `work`.
299 6334 : call setFFTF(factor, coef, g, f, inf)
300 6334 : if (inf) then
301 : ! fft of `g` is in `f`.
302 179858 : SET_CCF(f,work,f)
303 : else
304 : ! fft of `g` is in `g`.
305 0 : SET_CCF(f,work,g)
306 : end if
307 : else
308 : ! fft of `f` is in `f`.
309 5040 : call setFFTF(factor, coef, g, work, inf)
310 5040 : if (inf) then
311 : ! fft of `g` is in `work`.
312 0 : SET_CCF(f,f,work)
313 : else
314 : ! fft of `g` is in `g`.
315 143774 : SET_CCF(f,f,g)
316 : end if
317 : end if
318 : end block
319 : #else
320 : #error "Unrecognized interface."
321 : #endif
322 : ! result is now in `f`.
323 22744 : call setFFTR(factor, coef, f, g, inf)
324 22744 : inf = .not. inf
325 : !if (inwork) work = g
326 : !call setFFTR(factor, coef, work, f, inwork)
327 : !normfac = 1._TKC / real(size(work, 1, IK), TKC)
328 : !if (inwork) then
329 : ! do concurrent(isam = 1 : size(work, 1, IK))
330 : ! work(isam) = f(isam) * normfac
331 : ! end do
332 : !else
333 : ! do concurrent(isam = 1 : size(work, 1, IK))
334 : ! work(isam) = work(isam) * normfac
335 : ! end do
336 : !end if
337 : #undef SET_CCF
338 :
339 : #else
340 : !%%%%%%%%%%%%%%%%%%%%%%%%
341 : #error "Unrecognized interface."
342 : !%%%%%%%%%%%%%%%%%%%%%%%%
343 : #endif
344 : #undef TYPE_OF_SEQ
345 : #undef GET_ABSQ
346 : #undef GET_RE
|
