The ParaMonte Documentation Website
Current view: top level - kernel - CrossCorr_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: MPI Parallel Kernel - Code Coverage Report Lines: 309 325 95.1 %
Date: 2021-01-08 13:07:16 Functions: 16 16 100.0 %
Legend: Lines: hit not hit

          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        3615 :     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        2651 :         integer(IK)                         :: CumSumWeight(np), currentSampleEndLoc, batchSizeDefault, sampleCount, isample
      96             :         integer(IK)                         :: ip, ipVerbose, ipStart, ipEnd, npEffective
      97        2651 :         real(RK)                            :: avgPoint, varPoint, avgBatchMean, varBatchMean
      98        2651 :         real(RK)                            :: diffSquared, batchSizeDefaultInverse
      99             :         real(RK), allocatable               :: BatchMean(:)
     100             : 
     101             :         !Err%occurred = .false.
     102             : 
     103        2651 :         if (present(Weight)) then
     104         964 :             CumSumWeight = getCumSum(np,Weight)
     105             :         else
     106        1687 :             CumSumWeight(np) = np
     107             :         end if
     108             : 
     109             :         ! compute batch size and count
     110        2651 :         if (present(batchSize)) then
     111           3 :             batchSizeDefault = batchSize
     112             :         else
     113        2648 :             batchSizeDefault = int( real(CumSumWeight(np),kind=RK)**(0.666666666666666_RK) )
     114             :         end if
     115        2651 :         batchSizeDefaultInverse = 1._RK / real(batchSizeDefault,kind=RK)
     116        2651 :         sampleCount = CumSumWeight(np) / batchSizeDefault
     117        2651 :         npEffective = sampleCount * batchSizeDefault
     118             : 
     119        2651 :         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        2651 :         if (allocated(BatchMean)) deallocate(BatchMean)
     129        2651 :         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        2651 :         avgPoint = 0._RK
     136        2651 :         ipStart = 1
     137        2651 :         ipEnd   = 0
     138        2651 :         if (present(Weight)) then
     139         964 :             ip = 1
     140         964 :             isample = 1
     141         964 :             ipVerbose = 0
     142         964 :             currentSampleEndLoc = batchSizeDefault
     143         964 :             BatchMean(isample) = 0._RK
     144      491272 :             loopOverWeight: do
     145      492236 :                 ipVerbose = ipVerbose + 1
     146      492236 :                 if (ipVerbose>CumSumWeight(ip)) ip = ip + 1
     147      492236 :                 if (ipVerbose>currentSampleEndLoc) then ! we are done with the current batch
     148        6373 :                     avgPoint = avgPoint + BatchMean(isample)
     149        6373 :                     BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
     150        6373 :                     if (ipVerbose>npEffective) exit loopOverWeight  ! condition equivalent to currentSampleEndLoc==npEffective
     151        5409 :                     currentSampleEndLoc = currentSampleEndLoc + batchSizeDefault
     152        5409 :                     isample = isample + 1
     153        5409 :                     BatchMean(isample) = 0._RK
     154             :                 end if
     155      491272 :                 BatchMean(isample)  = BatchMean(isample) + Point(ip)
     156             :             end do loopOverWeight
     157             :         else    ! there is no weight
     158       10863 :             do isample = 1, sampleCount
     159        9176 :                 BatchMean(isample) = 0._RK
     160        9176 :                 ipEnd = ipEnd + batchSizeDefault
     161      408354 :                 do ip = ipStart, ipEnd
     162      408354 :                     BatchMean(isample)  = BatchMean(isample) + Point(ip)
     163             :                 end do
     164        9176 :                 ipStart = ipEnd + 1_IK
     165        9176 :                 avgPoint = avgPoint + BatchMean(isample)
     166       10863 :                 BatchMean(isample) = BatchMean(isample) * batchSizeDefaultInverse
     167             :             end do
     168             :         end if
     169       18200 :         avgBatchMean = sum( BatchMean ) / real(sampleCount,kind=RK)
     170       18200 :         varBatchMean = sum( (BatchMean - avgBatchMean)**2 ) / real(sampleCount-1,kind=RK)
     171        2651 :         avgPoint = avgPoint / real(npEffective,kind=RK)
     172             : 
     173             :         ! compute the variance of Point
     174             : 
     175        2651 :         varPoint = 0._RK
     176        2651 :         if (present(Weight)) then
     177         964 :             ip = 1
     178         964 :             ipVerbose = 0
     179         964 :             diffSquared = ( Point(ip) - avgPoint )**2
     180      491272 :             loopComputeVarPoint: do
     181      492236 :                 ipVerbose = ipVerbose + 1
     182      492236 :                 if (ipVerbose>npEffective) exit loopComputeVarPoint
     183      491272 :                 if (ipVerbose>CumSumWeight(ip)) then
     184      186645 :                     ip = ip + 1 ! by definition, ip never become > np, otherwise it leads to disastrous errors
     185      186645 :                     diffSquared = ( Point(ip) - avgPoint )**2
     186             :                 end if
     187      491272 :                 varPoint = varPoint + diffSquared
     188             :             end do loopComputeVarPoint
     189             :         else
     190      400865 :             do ip = 1, npEffective
     191      400865 :                 varPoint = varPoint + ( Point(ip) - avgPoint )**2
     192             :             end do
     193             :         end if
     194        2651 :         varPoint = varPoint / real(npEffective-1,kind=RK)
     195             : 
     196             :         ! compute the IAC
     197             : 
     198        2651 :         iac = batchSizeDefault * varBatchMean / varPoint
     199             : 
     200        2651 :     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         387 :     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        2651 :         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      217722 :         real(RK)                            :: iac, meanPoint, normFac, NormedData(np), cutoff, significanceDefault
     227             :         real(RK)    , allocatable           :: AutoCorr(:)
     228             :         integer(IK)                         :: i, paddedLen, sumWeight, cutoffIndex
     229         258 :         significanceDefault = 2._RK
     230         258 :         if (present(significance)) significanceDefault = significance
     231         258 :         if (present(Weight)) then
     232      117021 :             sumWeight = sum(Weight)
     233      117021 :             meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
     234             :         else
     235         129 :             sumWeight = np
     236      100701 :             meanPoint = sum(Point) / real(np,kind=RK)
     237             :         end if
     238      217722 :         NormedData = Point - meanPoint
     239         258 :         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         439 :                                             )
     248         258 :         normFac = 1._RK / AutoCorr(1)
     249     1623810 :         AutoCorr = AutoCorr * normFac
     250             :         ! For autocorrelation, under the assumption of a completely random series, the ACF standard error reduces to sqrt(1/ndata)
     251         258 :         cutoff = significanceDefault * sqrt(1._RK/sumWeight) ! standardErrorAutoCorr
     252         258 :         cutoffIndex = 1_IK
     253        7017 :         do i = 1, paddedLen
     254        7017 :             if (AutoCorr(i)<cutoff) then
     255         258 :                 cutoffIndex = i
     256         258 :                 exit
     257             :             end if
     258             :         end do
     259        7275 :         iac = 2_IK * sum(AutoCorr(1:cutoffIndex)) - 1_IK
     260         258 :     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         444 :     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         258 :         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      184053 :         real(RK)                            :: maxIAC, meanPoint, normFac, NormedData(np)
     284             :         real(RK), allocatable               :: AutoCorr(:)
     285             :         integer(IK)                         :: paddedLen, sumWeight
     286         303 :         if (present(Weight)) then
     287       78447 :             sumWeight = sum(Weight)
     288       78447 :             meanPoint = sum(Point*Weight) / real(sumWeight,kind=RK)
     289             :         else
     290         162 :             sumWeight = np
     291      105606 :             meanPoint = sum(Point) / real(np,kind=RK)
     292             :         end if
     293      184053 :         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         516 :                                             )
     303         303 :         normFac = 1._RK / AutoCorr(1)
     304     1174580 :         AutoCorr = AutoCorr * normFac
     305     1174880 :         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         576 :     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         576 :         baseDefault = 2_IK
     339           3 :         if (present(base)) baseDefault = base
     340         576 :         paddedLen = baseDefault ** ( getNextExponent( real(actualLen,kind=RK), real(baseDefault,kind=RK)) + 1_IK )
     341         576 :     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         576 :         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         600 :     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         600 :         if (present(base)) then
     428         588 :             nextExponent = ceiling( log(absoluteValue) / log(base) )
     429             :         else    ! assume the base is 2
     430          12 :             nextExponent = ceiling( log(absoluteValue) * INVLN2 )
     431             :         end if
     432         600 :     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         600 :         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     3292370 :     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     3291250 :         complex(CK), dimension(paddedLen/2)     :: Cdat1, Cdat2
     567             :         integer(IK)                             :: paddedLenHalf, paddedLenQuarter
     568             : 
     569         570 :         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         570 :         paddedLenHalf = paddedLen / 2
     575         570 :         paddedLenQuarter = paddedLenHalf / 2
     576         570 :         call realftWeighted(lenCompactData1, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData1, Cdat1, Weight1)
     577         570 :         call realftWeighted(lenCompactData2, paddedLen, paddedLenHalf, paddedLenQuarter, CompactData2, Cdat2, Weight2)
     578         570 :         Cdat1(1) = cmplx( real(Cdat1(1)) * real(Cdat2(1)) / paddedLenHalf, aimag(Cdat1(1)) * aimag(Cdat2(1)) / paddedLenHalf , kind = CK )
     579     1645060 :         Cdat1(2:paddedLenHalf) = Cdat1(2:paddedLenHalf) * conjg(Cdat2(2:paddedLenHalf)) / paddedLenHalf
     580         570 :         call realft(paddedLen, paddedLenHalf, paddedLenQuarter, CrossCorrFFT, -1_IK, Cdat1)
     581             : 
     582        1140 :     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        1158 :     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         570 :         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     1792510 :         complex(CK), dimension(paddedLenQuarter-1)  :: h1, h2
     601             :         complex(CK), pointer                        :: Cdata(:)
     602      896835 :         complex(CK)                                 :: w(paddedLenQuarter)
     603             :         complex(CK)                                 :: z
     604         579 :         real(RK)                                    :: c1, c2
     605             : 
     606         579 :         c1 = 0.5_RK
     607             : 
     608         579 :         if (present(Zdata)) then
     609         579 :             Cdata => Zdata
     610       98883 :             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         579 :         if (isign == 1_IK) then
     617           6 :             c2 = -0.5_RK
     618           6 :             call four1(paddedLenHalf, Cdata, +1_IK)
     619             :         else
     620         573 :             c2 = 0.5_RK
     621             :         end if
     622             : 
     623         579 :         w = zroots_unity(sign(paddedLen,isign), paddedLenQuarter)
     624      896835 :         w = cmplx(-aimag(w), real(w), kind=CK)
     625      896256 :         h1 = c1 * ( Cdata(2:paddedLenQuarter) + conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
     626      896256 :         h2 = c2 * ( Cdata(2:paddedLenQuarter) - conjg(Cdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
     627      896256 :         Cdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
     628      896256 :         Cdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
     629         579 :         z = Cdata(1)
     630             : 
     631         579 :         if (isign == 1_IK) then
     632           6 :             Cdata(1) = cmplx(real(z)+aimag(z), real(z)-aimag(z), kind=CK)
     633             :         else
     634         573 :             Cdata(1) = cmplx(c1*(real(z)+aimag(z)), c1*(real(z)-aimag(z)), kind=CK)
     635         573 :             call four1(paddedLenHalf, Cdata, -1)
     636             :         end if
     637             : 
     638         579 :      if (present(Zdata)) then
     639         579 :             if (isign /= 1_IK) then
     640     1694780 :                 PaddedData(1:paddedLen-1:2) = real(Cdata)
     641     1694780 :                 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         579 :     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        1692 :     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         579 :         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     1646200 :         complex(CK)                         :: w(paddedLenQuarter)
     670     1646200 :         complex(CK)                         :: h1(paddedLenQuarter-1)
     671     1646200 :         complex(CK)                         :: h2(paddedLenQuarter-1)
     672        1140 :         real(RK)                            :: zReal, zImag
     673        1140 :         real(RK)                            :: c1, c2
     674             : 
     675        1140 :         c1 = 0.5_RK
     676             : 
     677        1140 :         if (present(Weight)) then
     678             : 
     679         552 :             iw = 1
     680         552 :             counter = 0
     681      451176 :             loopOverCompactData: do id = 1, lenCompactData
     682      758166 :                 loopOverWeight: do
     683     1208800 :                     if (iw>Weight(id)) then
     684      225000 :                         iw = 1
     685      225000 :                         cycle loopOverCompactData
     686      983802 :                     elseif (iw==Weight(id)) then
     687      225906 :                         counter = counter + 1
     688      225906 :                         if (id==lenCompactData) then
     689         282 :                             Zdata(counter) = cmplx( CompactData(id) , 0._RK , kind=CK )
     690         282 :                             exit loopOverCompactData
     691             :                         else
     692      225624 :                             Zdata(counter) = cmplx( CompactData(id) , CompactData(id+1) , kind=CK )
     693      225624 :                             iw = 2
     694      225624 :                             cycle loopOverCompactData
     695             :                         end if
     696             :                     else
     697      757896 :                         counter = counter + 1
     698      757896 :                         Zdata(counter) = cmplx( CompactData(id) , CompactData(id) , kind=CK )
     699      757896 :                         iw = iw + 2
     700      757896 :                         cycle loopOverWeight
     701             :                     end if
     702             :                 end do loopOverWeight
     703             :             end do loopOverCompactData
     704     1566940 :             Zdata(counter+1:paddedLenHalf) = cmplx( 0._RK , 0._RK , kind=CK )
     705             : 
     706             :         else
     707             : 
     708         588 :             idEnd = lenCompactData / 2_IK
     709         588 :             do concurrent(id=1:idEnd)
     710      236490 :                 Zdata(id) = cmplx( CompactData(2*id-1) , CompactData(2*id) , kind=CK )
     711             :             end do
     712         588 :             if (2*idEnd<lenCompactData) then
     713         138 :                 idEnd = idEnd + 1
     714         138 :                 Zdata(idEnd) = cmplx( CompactData(lenCompactData) , 0._RK , kind=CK )
     715             :             end if
     716      504468 :             Zdata(idEnd+1:paddedLenHalf) = cmplx(0._RK, 0._RK, kind=CK)
     717             : 
     718             :         end if
     719             : 
     720        1140 :         c2 = -0.5_RK
     721        1140 :         call four1(paddedLenHalf,Zdata,1_IK)
     722             : 
     723        1140 :         w = zroots_unity(paddedLen, paddedLenQuarter)
     724     1646200 :         w = cmplx(-aimag(w), real(w,kind=RK), kind=CK)
     725     1645060 :         h1 = c1 * ( Zdata(2:paddedLenQuarter) + conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
     726     1645060 :         h2 = c2 * ( Zdata(2:paddedLenQuarter) - conjg(Zdata(paddedLenHalf:paddedLenQuarter+2:-1)) )
     727     1645060 :         Zdata(2:paddedLenQuarter) = h1 + w(2:paddedLenQuarter) * h2
     728     1645060 :         Zdata(paddedLenHalf:paddedLenQuarter+2:-1) = conjg(h1 - w(2:paddedLenQuarter) * h2)
     729        1140 :         zReal = real(Zdata(1)); zImag = aimag(Zdata(1))
     730        1140 :         Zdata(1) = cmplx(zReal+zImag, zReal-zImag, kind=CK)
     731             : 
     732        2280 :     end subroutine realftWeighted
     733             : 
     734             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     735             : 
     736        1719 :     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        1140 :         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        1719 :         m1=2**ceiling(0.5_RK*log(real(ndata,RK))/0.693147_RK)
     750        1719 :         m2=ndata/m1
     751        1719 :         allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
     752        5157 :         dat=reshape(Data,shape(dat))
     753        1719 :         call fourrow(dat,isign)
     754      127494 :         theta=arth(0,isign,m1)*TWOPI/ndata
     755      127494 :         wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
     756      125775 :         w=cmplx(1.0_RK,0.0_RK,kind=CK)
     757       34164 :         do j=2,m2
     758     5023460 :             w=w*wp+w
     759     4992730 :             dat(:,j)=dat(:,j)*w
     760             :         end do
     761     5210120 :         temp=transpose(dat)
     762        1719 :         call fourrow(temp,isign)
     763        3438 :         Data=reshape(temp,shape(Data))
     764        1719 :         deallocate(dat,w,wp,theta,temp)
     765        1719 :     end subroutine four1
     766             : 
     767             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     768             : 
     769        6876 :     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        1719 :         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        3438 :         real(RK)                                   :: theta
     780      161658 :         complex(CK), dimension(size(Data,1))       :: temp
     781             :         complex(CK)                                :: w,wp
     782             :         complex(CK)                                :: ws
     783        3438 :         n=size(Data,2)
     784        3438 :         n2=n/2
     785        3438 :         j=n2
     786      154854 :         do i=1,n-2
     787     4537220 :             if (j > i) call swap(Data(:,j+1),Data(:,i+1))
     788      151416 :             m=n2
     789      138654 :             do
     790      290070 :                 if (m < 2 .or. j < m) exit
     791      138654 :                 j=j-m
     792      138654 :                 m=m/2
     793             :             end do
     794      154854 :             j=j+m
     795             :         end do
     796        3438 :         mmax=1
     797       16128 :         do
     798       19566 :             if (n <= mmax) exit
     799       16128 :             istep=2*mmax
     800       16128 :             theta=PI/(isign*mmax)
     801       16128 :             wp=cmplx(-2.0_RK*sin(0.5_RK*theta)**2,sin(theta),kind=CK)
     802       16128 :             w=cmplx(1.0_RK,0.0_RK,kind=CK)
     803      170910 :             do m=1,mmax
     804      154782 :                 ws=w
     805      654228 :                 do i=m,n,istep
     806      499446 :                     j=i+mmax
     807    35311700 :                     temp=ws*Data(:,j)
     808    35311700 :                     Data(:,j)=Data(:,i)-temp
     809    35466500 :                     Data(:,i)=Data(:,i)+temp
     810             :                 end do
     811      170910 :                 w=w*wp+w
     812             :             end do
     813       16128 :             mmax=istep
     814             :         end do
     815        3438 :     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        3438 :         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

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021