!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a 
!!!!   copy of this software and associated documentation files (the "Software"), 
!!!!   to deal in the Software without restriction, including without limitation 
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense, 
!!!!   and/or sell copies of the Software, and to permit persons to whom the 
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be 
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR 
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE 
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of 
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your 
!!!!   work (education/research/industry/development/...) by citing the ParaMonte 
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> \brief This module contains procedures for computing various statistical correlation coefficients.
!> \author Amir Shahmoradi

module CorrCoef_mod

    use Constants_mod, only: RK, IK
    use Err_mod, only: Err_type
    implicit none

    character(*), parameter :: MODULE_NAME = "@CorrCoef_mod"

    type :: CorrCoefSpearman_type
        integer(IK)             :: ndata
        real(RK), allocatable   :: Data1(:), Data2(:)
        real(RK)                :: rho, rhoProb, dStarStar, dStarStarSignif, dStarStarProb
        type(Err_type)          :: Err
    contains
        procedure, nopass :: get => getCorrCoefSpearman
    end type CorrCoefSpearman_type

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Return the Spearman correlation coefficient of the two input data arrays.
    !> 
    !> @param[in] ndata             :   The length of the two input data arrays.
    !> @param[in] Data1             :   The input data array of size `(1:ndata)`.
    !> @param[in] Data2             :   The input data array of size `(1:ndata)`.
    !> @param[out] rho              :   The Spearman rank correlation coefficient.
    !> @param[out] rhoProb          :   The two-sided significance level of the deviation of `rho` from zero.
    !> @param[out] dStarStar        :   The sum-squared difference of ranks of the two data vectors.
    !> @param[out] dStarStarSignif  :   The number of standard deviations by which `dStarStar` deviates from its null-hypothesis expected value.
    !> @param[out] dStarStarProb    :   The two-sided significance level of the deviation represented by `dStarStarSignif`.
    !> @param[out] Err              :   An object of [Err_type](@ref err_mod::err_type) class, indicating whether an error has occurred.
    !> 
    !> \remark
    !> A small value of either `dStarStarProb` or `rhoProb` indicates a significant correlation.
    subroutine getCorrCoefSpearman(ndata,Data1,Data2,rho,rhoProb,dStarStar,dStarStarSignif,dStarStarProb,Err)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: getCorrCoefSpearman
#endif
        use Constants_mod, only: IK, RK, SQRT2, SPR
        use Statistics_mod, only: getBetaCDF
        use Sort_mod, only: sortAscending

        implicit none

        integer(IK)     , intent(in)    :: ndata
        real(RK)        , intent(in)    :: Data1(ndata), Data2(ndata)
        real(RK)        , intent(out)   :: rho, rhoProb
        real(RK)        , intent(out)   :: dStarStar, dStarStarSignif, dStarStarProb
        type(Err_type)  , intent(out)   :: Err

        character(*)    , parameter     :: PROCEDURE_NAME = MODULE_NAME//"@getCorrCoefSpearman"
        real(RK)                        :: aved,df,en,en3n,fac,sf,sg,t,vard
        real(RK)                        :: WorkSpace1(ndata), WorkSpace2(ndata)

        WorkSpace1 = Data1
        WorkSpace2 = Data2
        call sortAscending(ndata,WorkSpace1,WorkSpace2,Err)
        if (Err%occurred) then
        ! LCOV_EXCL_START
            Err%msg = PROCEDURE_NAME//Err%msg
            return
        end if
        ! LCOV_EXCL_STOP
        call crank(ndata,WorkSpace1,sf)

        call sortAscending(ndata,WorkSpace2,WorkSpace1,Err)
        if (Err%occurred) then
        ! LCOV_EXCL_START
            Err%msg = PROCEDURE_NAME//Err%msg
            return
        end if
        ! LCOV_EXCL_STOP
        call crank(ndata,WorkSpace2,sg)

        WorkSpace1 = WorkSpace1 - WorkSpace2
        dStarStar = dot_product(WorkSpace1,WorkSpace1)
        en = ndata
        en3n = en**3-en
        aved = en3n/6.0_RK-(sf+sg)/12.0_RK
        fac = (1.0_RK-sf/en3n)*(1.0_RK-sg/en3n)
        vard = ((en-1.0_RK)*en**2*(en+1.0_RK)**2/36.0_RK)*fac
        dStarStarSignif = (dStarStar-aved)/sqrt(vard)
        dStarStarProb = erfc( real( abs(dStarStarSignif)/SQRT2 , kind=SPR ) )
        rho = (1.0_RK-(6.0_RK/en3n)*(dStarStar+(sf+sg)/12.0_RK))/sqrt(fac)
        fac = (1.0_RK+rho)*(1.0_RK-rho)
        if (fac > 0.0) then
            t = rho*sqrt((en-2.0_RK)/fac)
            df = en-2.0_RK
            rhoProb = getBetaCDF(0.5_RK*df,0.5_RK,df/(df+t**2))
        else
            rhoProb = 0.0
        end if

    contains

        subroutine crank(ndata,Array,s)
            use Constants_mod, only: IK, RK
            use Misc_mod, only : arth, copyArray
            implicit none
            integer(IK), intent(in) :: ndata
            real(RK), intent(out) :: s
            real(RK), intent(inout) :: Array(ndata)
            integer(IK) :: i,ndum,nties
            integer(IK) :: Tstart(ndata), Tend(ndata), Tie(ndata), Idx(ndata)
            Idx = arth(1,1,ndata)
            Tie = merge(1,0,Array == eoshift(Array,-1))
            Tie(1) = 0
            Array = Idx(:)
            if (all(Tie == 0)) then
                s = 0.0
                return
            end if
            call copyArray( pack(Idx,Tie<eoshift(Tie,1)), Tstart, nties, ndum )
            Tend(1:nties) = pack(Idx(:),Tie(:)>eoshift(Tie(:),1))
            do i = 1,nties
                Array(Tstart(i):Tend(i)) = (Tstart(i)+Tend(i))/2.0_RK
            end do
            Tend(1:nties) = Tend(1:nties)-Tstart(1:nties)+1
            s = sum(Tend(1:nties)**3-Tend(1:nties))
        end subroutine crank

    end subroutine getCorrCoefSpearman

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module CorrCoef_mod ! LCOV_EXCL_LINE