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 3746 : 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 2716 : integer(IK) :: CumSumWeight(np), currentSampleEndLoc, batchSizeDefault, sampleCount, isample
96 : integer(IK) :: ip, ipVerbose, ipStart, ipEnd, npEffective
97 2716 : real(RK) :: avgPoint, varPoint, avgBatchMean, varBatchMean
98 2716 : real(RK) :: diffSquared, batchSizeDefaultInverse
99 : real(RK), allocatable :: BatchMean(:)
100 :
101 : !Err%occurred = .false.
102 :
103 2716 : if (present(Weight)) then
104 1030 : CumSumWeight = getCumSum(np,Weight)
105 : else
106 1686 : CumSumWeight(np) = np
107 : end if
108 :
109 : ! compute batch size and count
110 2716 : if (present(batchSize)) then
111 3 : batchSizeDefault = batchSize
112 : else
113 2713 : batchSizeDefault = int( real(CumSumWeight(np),kind=RK)**(0.666666666666666_RK) )
114 : end if
115 2716 : batchSizeDefaultInverse = 1._RK / real(batchSizeDefault,kind=RK)
116 2716 : sampleCount = CumSumWeight(np) / batchSizeDefault
117 2716 : npEffective = sampleCount * batchSizeDefault
118 :
119 2716 : if (sampleCount<2) then
120 : !Err%occurred = .true.
121 : !Err%msg = PROCEDURE_NAME // ": sampleCount<10: " // num2str(sampleCount)
122 : !iac = -huge(iac)
123 13 : iac = 1
124 13 : 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 2703 : if (allocated(BatchMean)) deallocate(BatchMean)
129 2703 : 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 2703 : avgPoint = 0._RK
136 2703 : ipStart = 1
137 2703 : ipEnd = 0
138 2703 : if (present(Weight)) then
139 1022 : ip = 1
140 1022 : isample = 1
141 1022 : ipVerbose = 0
142 1022 : currentSampleEndLoc = batchSizeDefault
143 1022 : BatchMean(isample) = 0._RK
144 513035 : loopOverWeight: do
145 514057 : ipVerbose = ipVerbose + 1
146 514057 : if (ipVerbose>CumSumWeight(ip)) ip = ip + 1
147 514057 : if (ipVerbose>currentSampleEndLoc) then ! we are done with the current batch
148 6694 : avgPoint = avgPoint + BatchMean(isample)
149 6694 : BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
150 6694 : if (ipVerbose>npEffective) exit loopOverWeight ! condition equivalent to currentSampleEndLoc==npEffective
151 5672 : currentSampleEndLoc = currentSampleEndLoc + batchSizeDefault
152 5672 : isample = isample + 1
153 5672 : BatchMean(isample) = 0._RK
154 : end if
155 513035 : BatchMean(isample) = BatchMean(isample) + Point(ip)
156 : end do loopOverWeight
157 : else ! there is no weight
158 10691 : do isample = 1, sampleCount
159 9010 : BatchMean(isample) = 0._RK
160 9010 : ipEnd = ipEnd + batchSizeDefault
161 400395 : do ip = ipStart, ipEnd
162 400395 : BatchMean(isample) = BatchMean(isample) + Point(ip)
163 : end do
164 9010 : ipStart = ipEnd + 1_IK
165 9010 : avgPoint = avgPoint + BatchMean(isample)
166 10691 : BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
167 : end do
168 : end if
169 18407 : avgBatchMean = sum( BatchMean ) / real(sampleCount,kind=RK)
170 18407 : varBatchMean = sum( (BatchMean - avgBatchMean)**2 ) / real(sampleCount-1,kind=RK)
171 2703 : avgPoint = avgPoint / real(npEffective,kind=RK)
172 :
173 : ! compute the variance of Point
174 :
175 2703 : varPoint = 0._RK
176 2703 : if (present(Weight)) then
177 1022 : ip = 1
178 1022 : ipVerbose = 0
179 1022 : diffSquared = ( Point(ip) - avgPoint )**2
180 513035 : loopComputeVarPoint: do
181 514057 : ipVerbose = ipVerbose + 1
182 514057 : if (ipVerbose>npEffective) exit loopComputeVarPoint
183 513035 : if (ipVerbose>CumSumWeight(ip)) then
184 198522 : ip = ip + 1 ! by definition, ip never become > np, otherwise it leads to disastrous errors
185 198522 : diffSquared = ( Point(ip) - avgPoint )**2
186 : end if
187 513035 : varPoint = varPoint + diffSquared
188 : end do loopComputeVarPoint
189 : else
190 393066 : do ip = 1, npEffective
191 393066 : varPoint = varPoint + ( Point(ip) - avgPoint )**2
192 : end do
193 : end if
194 2703 : varPoint = varPoint / real(npEffective-1,kind=RK)
195 :
196 : ! compute the IAC
197 :
198 2703 : iac = batchSizeDefault * varBatchMean / varPoint
199 :
200 2716 : 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 390 : 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 2716 : 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 216204 : real(RK) :: iac, meanPoint, normFac, NormedData(np), cutoff, significanceDefault
227 : real(RK) , allocatable :: AutoCorr(:)
228 : integer(IK) :: i, paddedLen, sumWeight, cutoffIndex
229 264 : significanceDefault = 2._RK
230 264 : if (present(significance)) significanceDefault = significance
231 264 : if (present(Weight)) then
232 115179 : sumWeight = sum(Weight)
233 115179 : meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
234 : else
235 138 : sumWeight = np
236 101025 : meanPoint = sum(Point) / real(np,kind=RK)
237 : end if
238 216204 : NormedData = Point - meanPoint
239 264 : 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 447 : )
248 264 : normFac = 1._RK / AutoCorr(1)
249 1623820 : AutoCorr = AutoCorr * normFac
250 : ! For autocorrelation, under the assumption of a completely random series, the ACF standard error reduces to sqrt(1/ndata)
251 264 : cutoff = significanceDefault * sqrt(1._RK/sumWeight) ! standardErrorAutoCorr
252 264 : cutoffIndex = 1_IK
253 6939 : do i = 1, paddedLen
254 6939 : if (AutoCorr(i)<cutoff) then
255 264 : cutoffIndex = i
256 264 : exit
257 : end if
258 : end do
259 7203 : iac = 2_IK * sum(AutoCorr(1:cutoffIndex)) - 1_IK
260 264 : 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 459 : 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 264 : 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 181272 : real(RK) :: maxIAC, meanPoint, normFac, NormedData(np)
284 : real(RK), allocatable :: AutoCorr(:)
285 : integer(IK) :: paddedLen, sumWeight
286 303 : if (present(Weight)) then
287 82173 : sumWeight = sum(Weight)
288 82173 : meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
289 : else
290 147 : sumWeight = np
291 99099 : meanPoint = sum(Point) / real(np,kind=RK)
292 : end if
293 181272 : NormedData = Point - meanPoint
294 303 : 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 513 : )
303 303 : normFac = 1._RK / AutoCorr(1)
304 1183580 : AutoCorr = AutoCorr * normFac
305 1183880 : maxIAC = 2_IK * maxval(getCumSum(paddedLen,AutoCorr)) - 1_IK
306 303 : 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 582 : 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 303 : 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 582 : baseDefault = 2_IK
339 3 : if (present(base)) baseDefault = base
340 582 : paddedLen = baseDefault ** ( getNextExponent( real(actualLen,kind=RK), real(baseDefault,kind=RK)) + 1_IK )
341 582 : 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 6 : 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 582 : use Constants_mod, only: IK, RK
369 : real(RK) , intent(in) :: actualLen
370 : real(RK) , intent(in), optional :: base
371 6 : real(RK) :: baseDefault
372 : integer(IK) :: paddedLen
373 6 : baseDefault = 2._RK
374 6 : if (present(base)) baseDefault = base
375 6 : paddedLen = nint( baseDefault ** ( getNextExponent(real(actualLen,kind=RK), baseDefault) + 1_IK ) )
376 6 : 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 12 : 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 6 : 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 12 : if (present(base)) then
400 3 : previousExponent = floor( log(absoluteValue) / log(base) )
401 : else ! assume the base is 2
402 9 : previousExponent = floor( log(absoluteValue) * INVLN2 )
403 : end if
404 12 : 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 606 : 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 12 : 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 606 : if (present(base)) then
428 594 : nextExponent = ceiling( log(absoluteValue) / log(base) )
429 : else ! assume the base is 2
430 12 : nextExponent = ceiling( log(absoluteValue) * INVLN2 )
431 : end if
432 606 : 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 9 : 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 606 : 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 9 : if (present(paddedLen)) then
460 6 : paddedSize = paddedLen
461 : else
462 3 : paddedSize = 2 ** ( getNextExponent(real(currentLen,kind=RK)) + 1 )
463 : end if
464 9 : allocate(ArrayPadded(paddedSize))
465 9 : do concurrent(i=1:currentLen)
466 29988 : ArrayPadded(i) = Array(i)
467 : end do
468 9 : do concurrent(i=currentLen+1:paddedSize)
469 68388 : ArrayPadded(i) = 0._RK
470 : end do
471 9 : 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 98313 : 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 9 : 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 98310 : complex(CK), dimension(paddedLen/2) :: Cdat1, Cdat2
511 : integer(IK) :: paddedLenHalf, paddedLenQuarter
512 :
513 3 : 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 3 : paddedLenHalf = paddedLen / 2
519 3 : paddedLenQuarter = paddedLenHalf / 2
520 3 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, PaddedData1, 1_IK, Cdat1)
521 3 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, PaddedData2, 1_IK, Cdat2)
522 3 : Cdat1(1) = cmplx( real(Cdat1(1)) * real(Cdat2(1)) / paddedLenHalf, aimag(Cdat1(1)) * aimag(Cdat2(1)) / paddedLenHalf , kind = CK )
523 49152 : Cdat1(2:paddedLenHalf) = Cdat1(2:paddedLenHalf) * conjg(Cdat2(2:paddedLenHalf)) / paddedLenHalf
524 3 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, CrossCorrFFT, -1_IK, Cdat1)
525 :
526 3 : 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 3301420 : 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 3 : 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 3300260 : complex(CK), dimension(paddedLen/2) :: Cdat1, Cdat2
567 : integer(IK) :: paddedLenHalf, paddedLenQuarter
568 :
569 576 : 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 576 : paddedLenHalf = paddedLen / 2
575 576 : paddedLenQuarter = paddedLenHalf / 2
576 576 : call realftWeighted(lenCompactData1, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData1, Cdat1, Weight1)
577 576 : call realftWeighted(lenCompactData2, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData2, Cdat2, Weight2)
578 576 : Cdat1(1) = cmplx( real(Cdat1(1)) * real(Cdat2(1)) / paddedLenHalf, aimag(Cdat1(1)) * aimag(Cdat2(1)) / paddedLenHalf , kind = CK )
579 1649560 : Cdat1(2:paddedLenHalf) = Cdat1(2:paddedLenHalf) * conjg(Cdat2(2:paddedLenHalf)) / paddedLenHalf
580 576 : call realft(paddedLen, paddedLenHalf, paddedLenQuarter, CrossCorrFFT, -1_IK, Cdat1)
581 :
582 1152 : 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 1170 : 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 576 : 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 1797010 : complex(CK), dimension(paddedLenQuarter-1) :: h1, h2
601 : complex(CK), pointer :: Cdata(:)
602 899091 : complex(CK) :: w(paddedLenQuarter)
603 : complex(CK) :: z
604 585 : real(RK) :: c1, c2
605 :
606 585 : c1 = 0.5_RK
607 :
608 585 : if (present(Zdata)) then
609 585 : Cdata => Zdata
610 98889 : 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 585 : if (isign == 1_IK) then
617 6 : c2 = -0.5_RK
618 6 : call four1(paddedLenHalf, Cdata, +1_IK)
619 : else
620 579 : c2 = 0.5_RK
621 : end if
622 :
623 585 : w = zroots_unity(sign(paddedLen,isign), paddedLenQuarter)
624 899091 : w = cmplx(-aimag(w), real(w), kind=CK)
625 898506 : h1 = c1 * ( Cdata(2:paddedLenQuarter) + conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
626 898506 : h2 = c2 * ( Cdata(2:paddedLenQuarter) - conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
627 898506 : Cdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
628 898506 : Cdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
629 585 : z = Cdata(1)
630 :
631 585 : if (isign == 1_IK) then
632 6 : Cdata(1) = cmplx(real(z)+aimag(z), real(z)-aimag(z), kind=CK)
633 : else
634 579 : Cdata(1) = cmplx(c1*(real(z)+aimag(z)), c1*(real(z)-aimag(z)), kind=CK)
635 579 : call four1(paddedLenHalf, Cdata, -1)
636 : end if
637 :
638 585 : if (present(Zdata)) then
639 585 : if (isign /= 1_IK) then
640 1699290 : PaddedData(1:paddedLen-1:2) = real(Cdata)
641 1699290 : 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 585 : 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 1728 : 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 585 : 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 1650710 : complex(CK) :: w(paddedLenQuarter)
670 1650710 : complex(CK) :: h1(paddedLenQuarter-1)
671 1650710 : complex(CK) :: h2(paddedLenQuarter-1)
672 1152 : real(RK) :: zReal, zImag
673 1152 : real(RK) :: c1, c2
674 :
675 1152 : c1 = 0.5_RK
676 :
677 1152 : if (present(Weight)) then
678 :
679 576 : iw = 1
680 576 : counter = 0
681 455046 : loopOverCompactData: do id = 1, lenCompactData
682 777810 : loopOverWeight: do
683 1232060 : if (iw>Weight(id)) then
684 226374 : iw = 1
685 226374 : cycle loopOverCompactData
686 1005690 : elseif (iw==Weight(id)) then
687 228276 : counter = counter + 1
688 228276 : if (id==lenCompactData) then
689 180 : Zdata(counter) = cmplx( CompactData(id) , 0._RK , kind=CK )
690 180 : exit loopOverCompactData
691 : else
692 228096 : Zdata(counter) = cmplx( CompactData(id) , CompactData(id+1) , kind=CK )
693 228096 : iw = 2
694 228096 : cycle loopOverCompactData
695 : end if
696 : else
697 777414 : counter = counter + 1
698 777414 : Zdata(counter) = cmplx( CompactData(id) , CompactData(id) , kind=CK )
699 777414 : iw = iw + 2
700 777414 : cycle loopOverWeight
701 : end if
702 : end do loopOverWeight
703 : end do loopOverCompactData
704 1583450 : Zdata(counter+1:paddedLenHalf) = cmplx( 0._RK , 0._RK , kind=CK )
705 :
706 : else
707 :
708 576 : idEnd = lenCompactData / 2_IK
709 576 : do concurrent(id=1:idEnd)
710 230304 : Zdata(id) = cmplx( CompactData(2*id-1) , CompactData(2*id) , kind=CK )
711 : end do
712 576 : if (2*idEnd<lenCompactData) then
713 132 : idEnd = idEnd + 1
714 132 : Zdata(idEnd) = cmplx( CompactData(lenCompactData) , 0._RK , kind=CK )
715 : end if
716 481260 : Zdata(idEnd+1:paddedLenHalf) = cmplx(0._RK, 0._RK, kind=CK)
717 :
718 : end if
719 :
720 1152 : c2 = -0.5_RK
721 1152 : call four1(paddedLenHalf,Zdata,1_IK)
722 :
723 1152 : w = zroots_unity(paddedLen, paddedLenQuarter)
724 1650710 : w = cmplx(-aimag(w), real(w,kind=RK), kind=CK)
725 1649560 : h1 = c1 * ( Zdata(2:paddedLenQuarter) + conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
726 1649560 : h2 = c2 * ( Zdata(2:paddedLenQuarter) - conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
727 1649560 : Zdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
728 1649560 : Zdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
729 1152 : zReal = real(Zdata(1)); zImag = aimag(Zdata(1))
730 1152 : Zdata(1) = cmplx(zReal+zImag, zReal-zImag, kind=CK)
731 :
732 2304 : end subroutine realftWeighted
733 :
734 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
735 :
736 1737 : 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 1152 : 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 1737 : m1=2**ceiling(0.5_RK*log(real(ndata,RK))/0.693147_RK)
750 1737 : m2=ndata/m1
751 1737 : allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
752 5211 : dat=reshape(Data,shape(dat))
753 1737 : call fourrow(dat,isign)
754 127710 : theta=arth(0,isign,m1)*TWOPI/ndata
755 127710 : wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
756 125973 : w=cmplx(1.0_RK,0.0_RK,kind=CK)
757 34155 : do j=2,m2
758 5036720 : w=w*wp+w
759 5006040 : dat(:,j)=dat(:,j)*w
760 : end do
761 5223830 : temp=transpose(dat)
762 1737 : call fourrow(temp,isign)
763 3474 : Data=reshape(temp,shape(Data))
764 1737 : deallocate(dat,w,wp,theta,temp)
765 1737 : end subroutine four1
766 :
767 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
768 :
769 6948 : 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 1737 : 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 3474 : real(RK) :: theta
780 161865 : complex(CK), dimension(size(Data,1)) :: temp
781 : complex(CK) :: w,wp
782 : complex(CK) :: ws
783 3474 : n=size(Data,2)
784 3474 : n2=n/2
785 3474 : j=n2
786 154980 : do i=1,n-2
787 4547020 : if (j > i) call swap(Data(:,j+1),Data(:,i+1))
788 151506 : m=n2
789 138654 : do
790 290160 : if (m < 2 .or. j < m) exit
791 138654 : j=j-m
792 138654 : m=m/2
793 : end do
794 154980 : j=j+m
795 : end do
796 3474 : mmax=1
797 16263 : do
798 19737 : if (n <= mmax) exit
799 16263 : istep=2*mmax
800 16263 : theta=PI/(isign*mmax)
801 16263 : wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
802 16263 : w=cmplx(1.0_RK,0.0_RK,kind=CK)
803 171180 : do m=1,mmax
804 154917 : ws=w
805 655119 : do i=m,n,istep
806 500202 : j=i+mmax
807 35431800 : temp=ws*Data(:,j)
808 35431800 : Data(:,j)=Data(:,i)-temp
809 35586700 : Data(:,i)=Data(:,i)+temp
810 : end do
811 171180 : w=w*wp+w
812 : end do
813 16263 : mmax=istep
814 : end do
815 3474 : end subroutine fourrow
816 :
817 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
818 : ! The following are for the AutoCorrelation computation using the conventional definition of correlation.
819 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
820 :
821 36 : 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 3474 : 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 18 : InverseSumNormedDataSq = 0._RK
832 89874 : do ip = 1, np
833 179739 : InverseSumNormedDataSq = InverseSumNormedDataSq + NormedData(1:nd,ip)**2
834 : end do
835 18 : InverseSumNormedDataSq = 1._RK / InverseSumNormedDataSq
836 9 : 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 6 : 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 9 : 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 12 : real(RK) :: InverseSumNormedDataSqDefault(nd)
871 : integer(IK) :: ip, ilag
872 :
873 96 : if (any(Lag>np-1_IK)) then
874 0 : AutoCorr = -huge(1._RK)
875 0 : return
876 : end if
877 :
878 6 : if (present(InverseSumNormedDataSq)) then
879 0 : InverseSumNormedDataSqDefault = InverseSumNormedDataSq
880 : else
881 6 : InverseSumNormedDataSqDefault = getInverseSumNormedDataSq(nd, np, NormedData)
882 : end if
883 :
884 : ! Now compute the non-normalized covariances
885 :
886 96 : do ilag = 1, nlag
887 96 : if (Lag(ilag)==0_IK) then
888 12 : AutoCorr(1:nd,ilag) = 1._RK
889 : else
890 168 : AutoCorr(1:nd,ilag) = 0._RK
891 740526 : do ip = 1, np-Lag(ilag)
892 1480970 : AutoCorr(1:nd,ilag) = AutoCorr(1:nd,ilag) + NormedData(1:nd,ip) * NormedData(1:nd,ip+Lag(ilag))
893 : end do
894 168 : AutoCorr(1:nd,ilag) = AutoCorr(1:nd,ilag) * InverseSumNormedDataSqDefault
895 : end if
896 : end do
897 :
898 12 : end subroutine getAutoCorrDirect
899 :
900 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
901 :
902 : end module CrossCorr_mod
|
