Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!!
4 : !!!! MIT License
5 : !!!!
6 : !!!! ParaMonte: plain powerful parallel Monte Carlo library.
7 : !!!!
8 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab
9 : !!!!
10 : !!!! This file is part of the ParaMonte library.
11 : !!!!
12 : !!!! Permission is hereby granted, free of charge, to any person obtaining a
13 : !!!! copy of this software and associated documentation files (the "Software"),
14 : !!!! to deal in the Software without restriction, including without limitation
15 : !!!! the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : !!!! and/or sell copies of the Software, and to permit persons to whom the
17 : !!!! Software is furnished to do so, subject to the following conditions:
18 : !!!!
19 : !!!! The above copyright notice and this permission notice shall be
20 : !!!! included in all copies or substantial portions of the Software.
21 : !!!!
22 : !!!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 : !!!! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 : !!!! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
25 : !!!! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
26 : !!!! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
27 : !!!! OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
28 : !!!! OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
29 : !!!!
30 : !!!! ACKNOWLEDGMENT
31 : !!!!
32 : !!!! ParaMonte is an honor-ware and its currency is acknowledgment and citations.
33 : !!!! As per the ParaMonte library license agreement terms, if you use any parts of
34 : !!!! this library for any purposes, kindly acknowledge the use of ParaMonte in your
35 : !!!! work (education/research/industry/development/...) by citing the ParaMonte
36 : !!!! library as described on this page:
37 : !!!!
38 : !!!! https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
39 : !!!!
40 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
42 :
43 : !> \brief This module contains procedures for computing the cross-correlation of time series data.
44 : !> \author Amir Shahmoradi
45 :
46 : module CrossCorr_mod
47 :
48 : implicit none
49 :
50 : character(len=*), parameter :: MODULE_NAME = "@CrossCorr_mod"
51 :
52 : interface getPaddedLen
53 : module procedure :: getPaddedLen_IK, getPaddedLen_RK
54 : end interface getPaddedLen
55 :
56 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 :
58 : contains
59 :
60 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 :
62 : !> \brief
63 : !> Return the integrated autocorrelation (IAC) via the BatchMeans method.
64 : !>
65 : !> @param[in] np : The number of data points in the input time series data.
66 : !> @param[in] Point : The input data series data vector.
67 : !> @param[in] Weight : The vector of weights of the input data points (**optional**, default = array of ones).
68 : !> @param[in] batchSize : The batch size (**optional**, default = computed from the input parameters).
69 : !>
70 : !> \return
71 : !> `iac` : The integrated autocorrelation (IAC) via the BatchMeans method.
72 : !>
73 : !> \remark
74 : !> Note that np must be large enough to get a meaningful answer.
75 3233 : function getBatchMeansIAC(np,Point,Weight,batchSize) result(iac)
76 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
77 : !DEC$ ATTRIBUTES DLLEXPORT :: getBatchMeansIAC
78 : #endif
79 : use Constants_mod, only: IK, RK
80 : use Statistics_mod, only: getVariance
81 : use Math_mod, only: getCumSum
82 : !use Stirng_mod, only: num2str
83 : !use Err_mod, only: Err_type
84 :
85 : implicit none
86 :
87 : character(len=*), parameter :: PROCEDURE_NAME = MODULE_NAME // "@getBatchMeansIAC()"
88 :
89 : integer(IK), intent(in) :: np
90 : real(RK), intent(in) :: Point(np)
91 : integer(IK), intent(in), optional :: Weight(np), batchSize
92 : !type(Err_type), intent(out) :: Err
93 : real(RK) :: iac
94 :
95 2376 : integer(IK) :: CumSumWeight(np), currentSampleEndLoc, batchSizeDefault, sampleCount, isample
96 : integer(IK) :: ip, ipVerbose, ipStart, ipEnd, npEffective
97 2376 : real(RK) :: avgPoint, varPoint, avgBatchMean, varBatchMean
98 2376 : real(RK) :: diffSquared, batchSizeDefaultInverse
99 : real(RK), allocatable :: BatchMean(:)
100 :
101 : !Err%occurred = .false.
102 :
103 2376 : if (present(Weight)) then
104 857 : CumSumWeight = getCumSum(np,Weight)
105 : else
106 1519 : CumSumWeight(np) = np
107 : end if
108 :
109 : ! compute batch size and count
110 2376 : if (present(batchSize)) then
111 1 : batchSizeDefault = batchSize
112 : else
113 2375 : batchSizeDefault = int( real(CumSumWeight(np),kind=RK)**(0.666666666666666_RK) )
114 : end if
115 2376 : batchSizeDefaultInverse = 1._RK / real(batchSizeDefault,kind=RK)
116 2376 : sampleCount = CumSumWeight(np) / batchSizeDefault
117 2376 : npEffective = sampleCount * batchSizeDefault
118 :
119 2376 : if (sampleCount<2) then
120 : !Err%occurred = .true.
121 : !Err%msg = PROCEDURE_NAME // ": sampleCount<10: " // num2str(sampleCount)
122 : !iac = -huge(iac)
123 0 : iac = 1
124 0 : return
125 : end if
126 :
127 : ! xxx: here goes another GFortran 7.3 bug: EndOfLineLoc is assumed already allocated, despite the first appearance here.
128 2376 : if (allocated(BatchMean)) deallocate(BatchMean)
129 2376 : allocate(BatchMean(sampleCount))
130 :
131 : ! improvement: iterate from the end to the beginning of the chain to ignore initial points instead of the last points.
132 : ! this would be beneficial for MCMC samples
133 :
134 : ! compute the Batch-Means avergage and variance, also average of Point
135 2376 : avgPoint = 0._RK
136 2376 : ipStart = 1
137 2376 : ipEnd = 0
138 2376 : if (present(Weight)) then
139 857 : ip = 1
140 857 : isample = 1
141 857 : ipVerbose = 0
142 857 : currentSampleEndLoc = batchSizeDefault
143 857 : BatchMean(isample) = 0._RK
144 280591 : loopOverWeight: do
145 281448 : ipVerbose = ipVerbose + 1
146 281448 : if (ipVerbose>CumSumWeight(ip)) ip = ip + 1
147 281448 : if (ipVerbose>currentSampleEndLoc) then ! we are done with the current batch
148 4852 : avgPoint = avgPoint + BatchMean(isample)
149 4852 : BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
150 4852 : if (ipVerbose>npEffective) exit loopOverWeight ! condition equivalent to currentSampleEndLoc==npEffective
151 3995 : currentSampleEndLoc = currentSampleEndLoc + batchSizeDefault
152 3995 : isample = isample + 1
153 3995 : BatchMean(isample) = 0._RK
154 : end if
155 280591 : BatchMean(isample) = BatchMean(isample) + Point(ip)
156 : end do loopOverWeight
157 : else ! there is no weight
158 9697 : do isample = 1, sampleCount
159 8178 : BatchMean(isample) = 0._RK
160 8178 : ipEnd = ipEnd + batchSizeDefault
161 342546 : do ip = ipStart, ipEnd
162 342546 : BatchMean(isample) = BatchMean(isample) + Point(ip)
163 : end do
164 8178 : ipStart = ipEnd + 1_IK
165 8178 : avgPoint = avgPoint + BatchMean(isample)
166 9697 : BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
167 : end do
168 : end if
169 15406 : avgBatchMean = sum( BatchMean ) / real(sampleCount,kind=RK)
170 15406 : varBatchMean = sum( (BatchMean - avgBatchMean)**2 ) / real(sampleCount-1,kind=RK)
171 2376 : avgPoint = avgPoint / real(npEffective,kind=RK)
172 :
173 : ! compute the variance of Point
174 :
175 2376 : varPoint = 0._RK
176 2376 : if (present(Weight)) then
177 857 : ip = 1
178 857 : ipVerbose = 0
179 857 : diffSquared = ( Point(ip) - avgPoint )**2
180 280591 : loopComputeVarPoint: do
181 281448 : ipVerbose = ipVerbose + 1
182 281448 : if (ipVerbose>npEffective) exit loopComputeVarPoint
183 280591 : if (ipVerbose>CumSumWeight(ip)) then
184 133911 : ip = ip + 1 ! by definition, ip never become > np, otherwise it leads to disastrous errors
185 133911 : diffSquared = ( Point(ip) - avgPoint )**2
186 : end if
187 280591 : varPoint = varPoint + diffSquared
188 : end do loopComputeVarPoint
189 : else
190 335887 : do ip = 1, npEffective
191 335887 : varPoint = varPoint + ( Point(ip) - avgPoint )**2
192 : end do
193 : end if
194 2376 : varPoint = varPoint / real(npEffective-1,kind=RK)
195 :
196 : ! compute the IAC
197 :
198 2376 : iac = batchSizeDefault * varBatchMean / varPoint
199 :
200 2376 : end function getBatchMeansIAC
201 :
202 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
203 :
204 : !> \brief
205 : !> Return the integrated autocorrelation (IAC) based on the cumulative autocorrelation.
206 : !>
207 : !> @param[in] np : The number of data points in the input time series data.
208 : !> @param[in] Point : The input data series data vector.
209 : !> @param[in] Weight : The vector of weights of the input data points (**optional**, default = array of ones).
210 : !> @param[in] significance : The significance in units of standard deviation below which the autocorrelation is
211 : !> considered noise (**optional**, default = 2).
212 : !>
213 : !> \return
214 : !> `iac` : The integrated autocorrelation (IAC).
215 179 : function getCumSumIAC(np,Point,Weight,significance) result(iac)
216 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
217 : !DEC$ ATTRIBUTES DLLEXPORT :: getCumSumIAC
218 : #endif
219 2376 : use Constants_mod, only: IK, RK
220 : use Math_mod, only: getCumSum
221 : implicit none
222 : integer(IK) , intent(in) :: np
223 : real(RK) , intent(in) :: Point(np)
224 : integer(IK) , intent(in), optional :: Weight(np)
225 : real(RK) , intent(in), optional :: significance ! in units of sigma
226 79212 : real(RK) :: iac, meanPoint, normFac, NormedData(np), cutoff, significanceDefault
227 : real(RK) , allocatable :: AutoCorr(:)
228 : integer(IK) :: i, paddedLen, sumWeight, cutoffIndex
229 123 : significanceDefault = 2._RK
230 123 : if (present(significance)) significanceDefault = significance
231 123 : if (present(Weight)) then
232 39235 : sumWeight = sum(Weight)
233 39235 : meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
234 : else
235 67 : sumWeight = np
236 39977 : meanPoint = sum(Point) / real(np,kind=RK)
237 : end if
238 79212 : NormedData = Point - meanPoint
239 123 : paddedLen = getPaddedLen(sumWeight)
240 0 : AutoCorr = getCrossCorrWeightedFFT ( lenCompactData1 = np &
241 : , lenCompactData2 = np &
242 : , paddedLen = paddedLen &
243 : , CompactData1 = NormedData &
244 : , CompactData2 = NormedData &
245 : , Weight1 = Weight &
246 : , Weight2 = Weight &
247 212 : )
248 123 : normFac = 1._RK / AutoCorr(1)
249 546939 : AutoCorr = AutoCorr * normFac
250 : ! For autocorrelation, under the assumption of a completely random series, the ACF standard error reduces to sqrt(1/ndata)
251 123 : cutoff = significanceDefault * sqrt(1._RK/sumWeight) ! standardErrorAutoCorr
252 123 : cutoffIndex = 1_IK
253 2520 : do i = 1, paddedLen
254 2520 : if (AutoCorr(i)<cutoff) then
255 123 : cutoffIndex = i
256 123 : exit
257 : end if
258 : end do
259 2643 : iac = 2_IK * sum(AutoCorr(1:cutoffIndex)) - 1_IK
260 123 : end function getCumSumIAC
261 :
262 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
263 :
264 : !> \brief
265 : !> Return the integrated autocorrelation (IAC) based on the maximum cumulative autocorrelation.
266 : !>
267 : !> @param[in] np : The number of data points in the input time series data.
268 : !> @param[in] Point : The input data series data vector.
269 : !> @param[in] Weight : The vector of weights of the input data points (**optional**, default = array of ones).
270 : !>
271 : !> \return
272 : !> `maxIAC` : The integrated autocorrelation (IAC).
273 156 : function getMaxCumSumIAC(np,Point,Weight) result(maxIAC)
274 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
275 : !DEC$ ATTRIBUTES DLLEXPORT :: getMaxCumSumIAC
276 : #endif
277 123 : use Constants_mod, only: IK, RK
278 : use Math_mod, only: getCumSum
279 : implicit none
280 : integer(IK), intent(in) :: np
281 : real(RK), intent(in) :: Point(np)
282 : integer(IK), intent(in), optional :: Weight(np)
283 63805 : real(RK) :: maxIAC, meanPoint, normFac, NormedData(np)
284 : real(RK), allocatable :: AutoCorr(:)
285 : integer(IK) :: paddedLen, sumWeight
286 119 : if (present(Weight)) then
287 21266 : sumWeight = sum(Weight)
288 21266 : meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
289 : else
290 82 : sumWeight = np
291 42539 : meanPoint = sum(Point) / real(np,kind=RK)
292 : end if
293 63805 : NormedData = Point - meanPoint
294 119 : paddedLen = getPaddedLen(sumWeight)
295 0 : AutoCorr = getCrossCorrWeightedFFT ( lenCompactData1 = np &
296 : , lenCompactData2 = np &
297 : , paddedLen = paddedLen &
298 : , CompactData1 = NormedData &
299 : , CompactData2 = NormedData &
300 : , Weight1 = Weight &
301 : , Weight2 = Weight &
302 207 : )
303 119 : normFac = 1._RK / AutoCorr(1)
304 371655 : AutoCorr = AutoCorr * normFac
305 371774 : maxIAC = 2_IK * maxval(getCumSum(paddedLen,AutoCorr)) - 1_IK
306 119 : end function getMaxCumSumIAC
307 :
308 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
309 :
310 : !> \brief
311 : !> Return the smallest length of a vector that is a power of `base` and at least `base**2` times larger than the input length `actualLen`.
312 : !>
313 : !> @param[in] actualLen : The input vector length.
314 : !> @param[in] base : The integer-valued base of the exponentiation (**optional**, default = `2`).
315 : !>
316 : !> \return
317 : !> `paddedLen` : The minimum power-of-`base` length given `actualLen`.
318 : !>
319 : !> \remark
320 : !> This method is used to compute the cross-correlation. Usage:
321 : !> `actualLen = max( actualLen1, actualLen2 )`
322 : !>
323 : !> \warning
324 : !> The input values for `absoluteValue` and `base` must be both larger than 1.
325 : !>
326 : !> \remark
327 : !> For weighted-data cross-correlation computation, try,
328 : !> `actualLen = max( sum(Weight1(1:actualLen1)), sum(Weight2(1:actualLen2)) )`.
329 247 : pure function getPaddedLen_IK(actualLen,base) result(paddedLen)
330 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
331 : !DEC$ ATTRIBUTES DLLEXPORT :: getPaddedLen_IK
332 : #endif
333 119 : use Constants_mod, only: IK, RK
334 : integer(IK) , intent(in) :: actualLen
335 : integer(IK) , intent(in), optional :: base
336 : integer(IK) :: baseDefault
337 : integer(IK) :: paddedLen
338 247 : baseDefault = 2_IK
339 1 : if (present(base)) baseDefault = base
340 247 : paddedLen = baseDefault ** ( getNextExponent( real(actualLen,kind=RK), real(baseDefault,kind=RK)) + 1_IK )
341 247 : end function getPaddedLen_IK
342 :
343 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
344 :
345 : !> \brief
346 : !> Return the smallest length of a vector that is a power of `base` and at least `base**2` times larger than the input length `actualLen`.
347 : !>
348 : !> @param[in] actualLen : The real-valued input vector length.
349 : !> @param[in] base : The real-valued base of the exponentiation (**optional**, default = `2.`).
350 : !>
351 : !> \return
352 : !> `paddedLen` : The minimum power-of-`base` length given `actualLen`.
353 : !>
354 : !> \remark
355 : !> This method is used to compute the cross-correlation. Usage:
356 : !> `actualLen = max( actualLen1, actualLen2 )`
357 : !>
358 : !> \warning
359 : !> The input values for `absoluteValue` and `base` must be both larger than 1.
360 : !>
361 : !> \remark
362 : !> For weighted-data cross-correlation computation, try,
363 : !> `actualLen = max( sum(Weight1(1:actualLen1)), sum(Weight2(1:actualLen2)) )`.
364 2 : pure function getPaddedLen_RK(actualLen,base) result(paddedLen)
365 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
366 : !DEC$ ATTRIBUTES DLLEXPORT :: getPaddedLen_RK
367 : #endif
368 247 : use Constants_mod, only: IK, RK
369 : real(RK) , intent(in) :: actualLen
370 : real(RK) , intent(in), optional :: base
371 2 : real(RK) :: baseDefault
372 : integer(IK) :: paddedLen
373 2 : baseDefault = 2._RK
374 2 : if (present(base)) baseDefault = base
375 2 : paddedLen = nint( baseDefault ** ( getNextExponent(real(actualLen,kind=RK), baseDefault) + 1_IK ) )
376 2 : end function getPaddedLen_RK
377 :
378 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
379 :
380 : !> \brief
381 : !> Return the exponent that yields the largest real number smaller than **or equal to** the input number `absoluteValue`.
382 : !>
383 : !> @param[in] absoluteValue : The input real number.
384 : !> @param[in] base : The base of the exponentiation (**optional**, default = `2`).
385 : !>
386 : !> \return
387 : !> `previousExponent` : The output minimum integer exponent.
388 : !>
389 : !> \warning
390 : !> The input values for `absoluteValue` and `base` must be both larger than 1.
391 4 : pure function getPreviousExponent(absoluteValue,base) result(previousExponent)
392 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
393 : !DEC$ ATTRIBUTES DLLEXPORT :: getPreviousExponent
394 : #endif
395 2 : use Constants_mod, only: IK, RK, INVLN2
396 : real(RK), intent(in) :: absoluteValue
397 : real(RK), intent(in), optional :: base
398 : integer(IK) :: previousExponent
399 4 : if (present(base)) then
400 1 : previousExponent = floor( log(absoluteValue) / log(base) )
401 : else ! assume the base is 2
402 3 : previousExponent = floor( log(absoluteValue) * INVLN2 )
403 : end if
404 4 : end function getPreviousExponent
405 :
406 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
407 :
408 : !> \brief
409 : !> Return the exponent that yields the smallest real number larger than **or equal to** the input number `absoluteValue`.
410 : !>
411 : !> @param[in] absoluteValue : The input real number.
412 : !> @param[in] base : The base of the exponentiation (**optional**, default = `2`).
413 : !>
414 : !> \warning
415 : !> The input values for `absoluteValue` and `base` must be both larger than 1.
416 : !>
417 : !> \return
418 : !> `nextExponent` : The output minimum integer exponent.
419 255 : pure function getNextExponent(absoluteValue,base) result(nextExponent)
420 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
421 : !DEC$ ATTRIBUTES DLLEXPORT :: getNextExponent
422 : #endif
423 4 : use Constants_mod, only: IK, RK, INVLN2
424 : real(RK), intent(in) :: absoluteValue
425 : real(RK), intent(in), optional :: base
426 : integer(IK) :: nextExponent
427 255 : if (present(base)) then
428 251 : nextExponent = ceiling( log(absoluteValue) / log(base) )
429 : else ! assume the base is 2
430 4 : nextExponent = ceiling( log(absoluteValue) * INVLN2 )
431 : end if
432 255 : end function getNextExponent
433 :
434 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
435 :
436 : !> \brief
437 : !> Return an array that is extended and padded with zeros for the requested length `paddedLen`.
438 : !>
439 : !> @param[in] currentLen : The length of the input array.
440 : !> @param[in] Array : The array to be extended.
441 : !> @param[in] paddedLen : The requested new length of the array.
442 : !>
443 : !> \return
444 : !> `ArrayPadded` : The output extended array, padded with zeros.
445 : !>
446 : !> \warning
447 : !> The input variable `paddedLen` must be a power of two, such that \f$2^\texttt{paddedLen}\f$
448 : !> represents the smallest integer larger than both `ndata1` and `ndata2`.
449 3 : pure function padZero(currentLen,Array,paddedLen) result(ArrayPadded)
450 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
451 : !DEC$ ATTRIBUTES DLLEXPORT :: padZero
452 : #endif
453 255 : use Constants_mod, only: IK, RK
454 : integer(IK) , intent(in) :: currentLen
455 : real(RK) , intent(in) :: Array(currentLen)
456 : integer(IK) , intent(in), optional :: paddedLen
457 : real(RK) , allocatable :: ArrayPadded(:)
458 : integer(IK) :: i, paddedSize
459 3 : if (present(paddedLen)) then
460 2 : paddedSize = paddedLen
461 : else
462 1 : paddedSize = 2 ** ( getNextExponent(real(currentLen,kind=RK)) + 1 )
463 : end if
464 3 : allocate(ArrayPadded(paddedSize))
465 3 : do concurrent(i=1:currentLen)
466 9996 : ArrayPadded(i) = Array(i)
467 : end do
468 3 : do concurrent(i=currentLen+1:paddedSize)
469 22796 : ArrayPadded(i) = 0._RK
470 : end do
471 3 : end function padZero
472 :
473 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
474 :
475 : !> \brief
476 : !> Return the cross-correlation of the two input data vectors,
477 : !> (including any user-supplied zero padding), computed via Fast-Fourier Transform.
478 : !>
479 : !> @param[in] paddedLen : The lengths of the input arrays, which **MUST** be an integer *power of two*.
480 : !> @param[in] PaddedData1 : The first array, possibly padded with zero to become an array of length `paddedLen`.
481 : !> @param[in] PaddedData2 : The second array, possibly padded with zero to become an array of length `paddedLen`.
482 : !>
483 : !> \return
484 : !> `CrossCorrFFT` : A vector of same length as the input arrays containing the cross-correlation.
485 : !> The answer is returned as the first `paddedLen` points in ans stored in wrap-around order, i.e., cross-correlation
486 : !> at increasingly negative lags are in `CrossCorrFFT(paddedLen)` on down to `CrossCorrFFT(paddedLen/2+1)`, while
487 : !> cross-correlation increasingly positive lags are in `CrossCorrFFT(1)` (zero lag) on up to `CrossCorrFFT(paddedLen/2)`.
488 : !>
489 : !> \remark
490 : !> Sign convention of this routine: If `PaddedData1` lags `PaddedData2`, i.e.,
491 : !> is shifted to the right of it, then ans will show a peak at positive lags.
492 : !>
493 : !> \remark
494 : !> For autocorrelation, under the assumption of a completely random series,
495 : !> the ACF standard error reduces to: \f$\sqrt{ 1 / \texttt{paddedLen} }\f$
496 32771 : function getCrossCorrFFT(paddedLen, PaddedData1, PaddedData2) result(CrossCorrFFT)
497 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
498 : !DEC$ ATTRIBUTES DLLEXPORT :: getCrossCorrFFT
499 : #endif
500 : use iso_fortran_env, only: output_unit
501 3 : use Constants_mod, only: IK, RK, CK
502 :
503 : implicit none
504 :
505 : integer(IK), intent(in) :: paddedLen
506 : real(RK), intent(inout) :: PaddedData1(paddedLen), PaddedData2(paddedLen)
507 : real(RK) :: CrossCorrFFT(paddedLen)
508 :
509 : character(len=*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@getCrossCorrFFT()"
510 32770 : complex(CK), dimension(paddedLen/2) :: Cdat1, Cdat2
511 : integer(IK) :: paddedLenHalf, paddedLenQuarter
512 :
513 1 : if (iand(paddedLen,paddedLen-1) /= 0) then
514 0 : write(output_unit,"(A)") PROCEDURE_NAME//": paddedLen must be a power of 2."
515 0 : error stop
516 : end if
517 :
518 1 : paddedLenHalf = paddedLen / 2
519 1 : paddedLenQuarter = paddedLenHalf / 2
520 1 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, PaddedData1, 1_IK, Cdat1)
521 1 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, PaddedData2, 1_IK, Cdat2)
522 1 : Cdat1(1) = cmplx( real(Cdat1(1)) * real(Cdat2(1)) / paddedLenHalf, aimag(Cdat1(1)) * aimag(Cdat2(1)) / paddedLenHalf , kind = CK )
523 16384 : Cdat1(2:paddedLenHalf) = Cdat1(2:paddedLenHalf) * conjg(Cdat2(2:paddedLenHalf)) / paddedLenHalf
524 1 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, CrossCorrFFT, -1_IK, Cdat1)
525 :
526 1 : end function getCrossCorrFFT
527 :
528 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
529 :
530 : !> \brief
531 : !> Return the cross-correlation of the two input *weighted* compact data vectors,
532 : !> (including any user-supplied zero padding), computed via Fast-Fourier Transform.
533 : !>
534 : !> @param[in] lenCompactData1 : The length of the first input array `CompactData1`.
535 : !> @param[in] lenCompactData2 : The length of the second input array `CompactData2`.
536 : !> @param[in] paddedLen : The length by which the input data vectors must be extended and padded.
537 : !> @param[in] CompactData1 : The first array, in compact weighted format.
538 : !> @param[in] CompactData2 : The second array, in compact weighted format.
539 : !> @param[in] Weight1 : The weights of the elements in the first array.
540 : !> @param[in] Weight2 : The weights of the elements in the second array.
541 : !>
542 : !> \return
543 : !> `CrossCorrFFT` : A vector of length `paddedLen` containing the cross-correlation.
544 : !>
545 : !> \warning
546 : !> The input variable `paddedLen` must be a power of two, such that \f$2^\texttt{paddedLen}\f$
547 : !> represents the smallest integer larger than both `lenCompactData1` and `lenCompactData2`.
548 : !>
549 : !> \remark
550 : !> Unlike [getCrossCorrFFT](@ref getcrosscorrfft), the input arguments `CompactData1`
551 : !> and `CompactData2` to this function can be the same arrays, as the `intent` of both is `in`.
552 1083370 : function getCrossCorrWeightedFFT(lenCompactData1, lenCompactData2, paddedLen, CompactData1, CompactData2, Weight1, Weight2) result(CrossCorrFFT)
553 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
554 : !DEC$ ATTRIBUTES DLLEXPORT :: getCrossCorrWeightedFFT
555 : #endif
556 : use iso_fortran_env, only: output_unit
557 1 : use Constants_mod, only: IK, RK, CK
558 : implicit none
559 : integer(IK), intent(in) :: lenCompactData1, lenCompactData2, paddedLen
560 : real(RK), intent(in) :: CompactData1(lenCompactData1), CompactData2(lenCompactData2)
561 : integer(IK), intent(in), optional :: Weight1(lenCompactData1), Weight2(lenCompactData2)
562 : real(RK) :: CrossCorrFFT(paddedLen)
563 :
564 : character(len=*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@getCrossCorrWeightedFFT()"
565 :
566 1082940 : complex(CK), dimension(paddedLen/2) :: Cdat1, Cdat2
567 : integer(IK) :: paddedLenHalf, paddedLenQuarter
568 :
569 245 : if (iand(paddedLen,paddedLen-1) /= 0) then
570 0 : write(output_unit,"(A)") PROCEDURE_NAME//": paddedLen must be a power of 2."
571 0 : error stop
572 : end if
573 :
574 245 : paddedLenHalf = paddedLen / 2
575 245 : paddedLenQuarter = paddedLenHalf / 2
576 245 : call realftWeighted(lenCompactData1, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData1, Cdat1, Weight1)
577 245 : call realftWeighted(lenCompactData2, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData2, Cdat2, Weight2)
578 245 : Cdat1(1) = cmplx( real(Cdat1(1)) * real(Cdat2(1)) / paddedLenHalf, aimag(Cdat1(1)) * aimag(Cdat2(1)) / paddedLenHalf , kind = CK )
579 541224 : Cdat1(2:paddedLenHalf) = Cdat1(2:paddedLenHalf) * conjg(Cdat2(2:paddedLenHalf)) / paddedLenHalf
580 245 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, CrossCorrFFT, -1_IK, Cdat1)
581 :
582 490 : end function getCrossCorrWeightedFFT
583 :
584 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
585 :
586 : ! Here, `paddedLenQuarter` is 1/4 of the length of the input array, which is already padded with zeros.
587 496 : subroutine realft(paddedLen, paddedLenHalf, paddedLenQuarter, PaddedData, isign, Zdata)
588 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
589 : !DEC$ ATTRIBUTES DLLEXPORT :: realft
590 : #endif
591 245 : use Constants_mod, only: IK, RK, CK
592 : use Misc_mod, only: zroots_unity
593 :
594 : implicit none
595 :
596 : integer(IK), intent(in) :: paddedLen, paddedLenHalf, paddedLenQuarter, isign
597 : real(RK), intent(inout) :: PaddedData(paddedLen)
598 : complex(CK), optional, target :: Zdata(paddedLenHalf)
599 :
600 590376 : complex(CK), dimension(paddedLenQuarter-1) :: h1, h2
601 : complex(CK), pointer :: Cdata(:)
602 295436 : complex(CK) :: w(paddedLenQuarter)
603 : complex(CK) :: z
604 248 : real(RK) :: c1, c2
605 :
606 248 : c1 = 0.5_RK
607 :
608 248 : if (present(Zdata)) then
609 248 : Cdata => Zdata
610 33016 : if (isign == 1_IK) Cdata = cmplx(PaddedData(1:paddedLen-1:2), PaddedData(2:paddedLen:2), kind=CK)
611 : else
612 0 : allocate(Cdata(paddedLenHalf))
613 0 : Cdata = cmplx(PaddedData(1:paddedLen-1:2), PaddedData(2:paddedLen:2), kind=CK)
614 : end if
615 :
616 248 : if (isign == 1_IK) then
617 2 : c2 = -0.5_RK
618 2 : call four1(paddedLenHalf, Cdata, +1_IK)
619 : else
620 246 : c2 = 0.5_RK
621 : end if
622 :
623 248 : w = zroots_unity(sign(paddedLen,isign), paddedLenQuarter)
624 295436 : w = cmplx(-aimag(w), real(w), kind=CK)
625 295188 : h1 = c1 * ( Cdata(2:paddedLenQuarter) + conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
626 295188 : h2 = c2 * ( Cdata(2:paddedLenQuarter) - conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
627 295188 : Cdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
628 295188 : Cdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
629 248 : z = Cdata(1)
630 :
631 248 : if (isign == 1_IK) then
632 2 : Cdata(1) = cmplx(real(z)+aimag(z), real(z)-aimag(z), kind=CK)
633 : else
634 246 : Cdata(1) = cmplx(c1*(real(z)+aimag(z)), c1*(real(z)-aimag(z)), kind=CK)
635 246 : call four1(paddedLenHalf, Cdata, -1)
636 : end if
637 :
638 248 : if (present(Zdata)) then
639 248 : if (isign /= 1_IK) then
640 557854 : PaddedData(1:paddedLen-1:2) = real(Cdata)
641 557854 : PaddedData(2:paddedLen:2) = aimag(Cdata)
642 : end if
643 : else
644 0 : PaddedData(1:paddedLen-1:2) = real(Cdata)
645 0 : PaddedData(2:paddedLen:2) = aimag(Cdata)
646 0 : nullify(Cdata)
647 : end if
648 :
649 248 : end subroutine realft
650 :
651 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
652 :
653 : ! The `lenCompactData` is the length of the compact input array `CompactData`, which is NOT padded by default.
654 680 : subroutine realftWeighted(lenCompactData, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData, Zdata, Weight)
655 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
656 : !DEC$ ATTRIBUTES DLLEXPORT :: realftWeighted
657 : #endif
658 248 : use Constants_mod, only: IK, RK, CK
659 : use Misc_mod, only: zroots_unity
660 :
661 : implicit none
662 :
663 : integer(IK) , intent(in) :: lenCompactData, paddedLen, paddedLenHalf, paddedLenQuarter
664 : real(RK) , intent(in) :: CompactData(lenCompactData)
665 : integer(IK) , intent(in), optional :: Weight(lenCompactData)
666 : complex(CK) , intent(out) :: Zdata(paddedLenHalf)
667 :
668 : integer(IK) :: id, idEnd, iw, counter
669 541714 : complex(CK) :: w(paddedLenQuarter)
670 541714 : complex(CK) :: h1(paddedLenQuarter-1)
671 541714 : complex(CK) :: h2(paddedLenQuarter-1)
672 490 : real(RK) :: zReal, zImag
673 490 : real(RK) :: c1, c2
674 :
675 490 : c1 = 0.5_RK
676 :
677 490 : if (present(Weight)) then
678 :
679 190 : iw = 1
680 190 : counter = 0
681 141050 : loopOverCompactData: do id = 1, lenCompactData
682 232902 : loopOverWeight: do
683 373824 : if (iw>Weight(id)) then
684 69968 : iw = 1
685 69968 : cycle loopOverCompactData
686 303856 : elseif (iw==Weight(id)) then
687 71018 : counter = counter + 1
688 71018 : if (id==lenCompactData) then
689 126 : Zdata(counter) = cmplx( CompactData(id) , 0._RK , kind=CK )
690 126 : exit loopOverCompactData
691 : else
692 70892 : Zdata(counter) = cmplx( CompactData(id) , CompactData(id+1) , kind=CK )
693 70892 : iw = 2
694 70892 : cycle loopOverCompactData
695 : end if
696 : else
697 232838 : counter = counter + 1
698 232838 : Zdata(counter) = cmplx( CompactData(id) , CompactData(id) , kind=CK )
699 232838 : iw = iw + 2
700 232838 : cycle loopOverWeight
701 : end if
702 : end do loopOverWeight
703 : end do loopOverCompactData
704 486894 : Zdata(counter+1:paddedLenHalf) = cmplx( 0._RK , 0._RK , kind=CK )
705 :
706 : else
707 :
708 300 : idEnd = lenCompactData / 2_IK
709 300 : do concurrent(id=1:idEnd)
710 92628 : Zdata(id) = cmplx( CompactData(2*id-1) , CompactData(2*id) , kind=CK )
711 : end do
712 300 : if (2*idEnd<lenCompactData) then
713 48 : idEnd = idEnd + 1
714 48 : Zdata(idEnd) = cmplx( CompactData(lenCompactData) , 0._RK , kind=CK )
715 : end if
716 199812 : Zdata(idEnd+1:paddedLenHalf) = cmplx(0._RK, 0._RK, kind=CK)
717 :
718 : end if
719 :
720 490 : c2 = -0.5_RK
721 490 : call four1(paddedLenHalf,Zdata,1_IK)
722 :
723 490 : w = zroots_unity(paddedLen, paddedLenQuarter)
724 541714 : w = cmplx(-aimag(w), real(w,kind=RK), kind=CK)
725 541224 : h1 = c1 * ( Zdata(2:paddedLenQuarter) + conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
726 541224 : h2 = c2 * ( Zdata(2:paddedLenQuarter) - conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
727 541224 : Zdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
728 541224 : Zdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
729 490 : zReal = real(Zdata(1)); zImag = aimag(Zdata(1))
730 490 : Zdata(1) = cmplx(zReal+zImag, zReal-zImag, kind=CK)
731 :
732 980 : end subroutine realftWeighted
733 :
734 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
735 :
736 738 : subroutine four1(ndata,Data,isign)
737 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
738 : !DEC$ ATTRIBUTES DLLEXPORT :: four1
739 : #endif
740 490 : use Constants_mod, only: IK, RK, CK, TWOPI
741 : use Misc_mod, only: arth
742 : implicit none
743 : integer(IK), intent(in) :: ndata, isign
744 : complex(CK), intent(inout) :: Data(ndata)
745 : complex(CK), allocatable :: dat(:,:), temp(:,:)
746 : complex(CK), allocatable :: w(:),wp(:)
747 : real(RK) , allocatable :: theta(:)
748 : integer(IK) :: m1,m2,j
749 738 : m1=2**ceiling(0.5_RK*log(real(ndata,RK))/0.693147_RK)
750 738 : m2=ndata/m1
751 738 : allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
752 2214 : dat=reshape(Data,shape(dat))
753 738 : call fourrow(dat,isign)
754 46092 : theta=arth(0,isign,m1)*TWOPI/ndata
755 46092 : wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
756 45354 : w=cmplx(1.0_RK,0.0_RK,kind=CK)
757 12696 : do j=2,m2
758 1652120 : w=w*wp+w
759 1640900 : dat(:,j)=dat(:,j)*w
760 : end do
761 1718920 : temp=transpose(dat)
762 738 : call fourrow(temp,isign)
763 1476 : Data=reshape(temp,shape(Data))
764 738 : deallocate(dat,w,wp,theta,temp)
765 738 : end subroutine four1
766 :
767 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
768 :
769 2952 : subroutine fourrow(Data,isign)
770 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
771 : !DEC$ ATTRIBUTES DLLEXPORT :: fourrow
772 : #endif
773 738 : use Constants_mod, only: IK, RK, CK, PI
774 : use Misc_mod, only: swap
775 : implicit none
776 : complex(CK), dimension(:,:), intent(inout) :: Data
777 : integer(IK), intent(in) :: isign
778 : integer(IK) :: n,i,istep,j,m,mmax,n2
779 1476 : real(RK) :: theta
780 58788 : complex(CK), dimension(size(Data,1)) :: temp
781 : complex(CK) :: w,wp
782 : complex(CK) :: ws
783 1476 : n=size(Data,2)
784 1476 : n2=n/2
785 1476 : j=n2
786 55854 : do i=1,n-2
787 1492610 : if (j > i) call swap(Data(:,j+1),Data(:,i+1))
788 54378 : m=n2
789 49146 : do
790 103524 : if (m < 2 .or. j < m) exit
791 49146 : j=j-m
792 49146 : m=m/2
793 : end do
794 55854 : j=j+m
795 : end do
796 1476 : mmax=1
797 6690 : do
798 8166 : if (n <= mmax) exit
799 6690 : istep=2*mmax
800 6690 : theta=PI/(isign*mmax)
801 6690 : wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
802 6690 : w=cmplx(1.0_RK,0.0_RK,kind=CK)
803 62526 : do m=1,mmax
804 55836 : ws=w
805 229593 : do i=m,n,istep
806 173757 : j=i+mmax
807 11538400 : temp=ws*Data(:,j)
808 11538400 : Data(:,j)=Data(:,i)-temp
809 11594200 : Data(:,i)=Data(:,i)+temp
810 : end do
811 62526 : w=w*wp+w
812 : end do
813 6690 : mmax=istep
814 : end do
815 1476 : end subroutine fourrow
816 :
817 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
818 : ! The following are for the AutoCorrelation computation using the conventional definition of correlation.
819 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
820 :
821 12 : pure function getInverseSumNormedDataSq(nd,np,NormedData) result(InverseSumNormedDataSq)
822 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
823 : !DEC$ ATTRIBUTES DLLEXPORT :: getInverseSumNormedDataSq
824 : #endif
825 1476 : use Constants_mod, only: IK, RK
826 : implicit none
827 : integer(IK), intent(in) :: nd,np
828 : real(RK) , intent(in) :: NormedData(nd,np)
829 : real(RK) :: InverseSumNormedDataSq(nd)
830 : integer(IK) :: ip
831 6 : InverseSumNormedDataSq = 0._RK
832 29958 : do ip = 1, np
833 59913 : InverseSumNormedDataSq = InverseSumNormedDataSq + NormedData(1:nd,ip)**2
834 : end do
835 6 : InverseSumNormedDataSq = 1._RK / InverseSumNormedDataSq
836 3 : end function getInverseSumNormedDataSq
837 :
838 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
839 :
840 : !> \brief
841 : !> Compute the autocorrelation of the input data matrix (that is already normalized with respect to its mean).
842 : !>
843 : !> \param[in] nd : The number of dimensions of the input data (the number of data features).
844 : !> \param[in] np : The length of the input data (the number of observations or points).
845 : !> \param[in] NormedData : The input data of size `(nd,np)`. On input, it must be normalized
846 : !> with respect to it mean vector.
847 : !> \param[in] nlag : The length of the input vector of lags.
848 : !> \param[in] Lag : The input vector of lags at which the autocorrelation must be computed.
849 : !> \param[out] AutoCorr : The output AutoCorrelation matrix of shape `(nd,nlag)`.
850 : !> \param[in] InverseSumNormedDataSq : The inverse sum of the square of the input normalized data (**optional**).
851 : !>
852 : !> \remark
853 : !> If `InverseSumNormedDataSq` is provided, the computations will be slightly faster.
854 : !>
855 : !> \remark
856 : !> This function uses the direct method of covariance computation for computing the autocorrelation.
857 : !> It is therefore highly inefficient and not recommended for use in HPC solutions.
858 : !> Use instead, [getCrossCorrFFT](@ref getcrosscorrfft) or [getCrossCorrWeightedFFT](@ref getcrosscorrweightedfft).
859 2 : pure subroutine getAutoCorrDirect(nd, np, NormedData, nlag, Lag, AutoCorr, InverseSumNormedDataSq)
860 : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
861 : !DEC$ ATTRIBUTES DLLEXPORT :: getAutoCorrDirect
862 : #endif
863 :
864 3 : use Constants_mod, only: IK, RK
865 : implicit none
866 : integer(IK), intent(in) :: nd,np,nlag,Lag(nlag)
867 : real(RK) , intent(in) :: NormedData(nd,np)
868 : real(RK) , intent(in), optional :: InverseSumNormedDataSq(nd)
869 : real(RK) , intent(out) :: AutoCorr(nd,nlag)
870 4 : real(RK) :: InverseSumNormedDataSqDefault(nd)
871 : integer(IK) :: ip, ilag
872 :
873 32 : if (any(Lag>np-1_IK)) then
874 0 : AutoCorr = -huge(1._RK)
875 0 : return
876 : end if
877 :
878 2 : if (present(InverseSumNormedDataSq)) then
879 0 : InverseSumNormedDataSqDefault = InverseSumNormedDataSq
880 : else
881 2 : InverseSumNormedDataSqDefault = getInverseSumNormedDataSq(nd, np, NormedData)
882 : end if
883 :
884 : ! Now compute the non-normalized covariances
885 :
886 32 : do ilag = 1, nlag
887 32 : if (Lag(ilag)==0_IK) then
888 4 : AutoCorr(1:nd,ilag) = 1._RK
889 : else
890 56 : AutoCorr(1:nd,ilag) = 0._RK
891 246842 : do ip = 1, np-Lag(ilag)
892 493656 : AutoCorr(1:nd,ilag) = AutoCorr(1:nd,ilag) + NormedData(1:nd,ip) * NormedData(1:nd,ip+Lag(ilag))
893 : end do
894 56 : AutoCorr(1:nd,ilag) = AutoCorr(1:nd,ilag) * InverseSumNormedDataSqDefault
895 : end if
896 : end do
897 :
898 4 : end subroutine getAutoCorrDirect
899 :
900 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
901 :
902 : end module CrossCorr_mod
|
