#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_stats) stdlib_stats_corr
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_error, only: error_stop
use stdlib_linalg, only: diag
use stdlib_optval, only: optval
implicit none
contains
#:for k1, t1 in RC_KINDS_TYPES
#:set RName = rname("corr",1, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:)
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(${k1}$) :: res
if (.not.optval(mask, .true.) .or. size(x) < 2) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if
res = 1
end function ${RName}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:set RName = rname("corr",1, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:)
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(dp) :: res
if (.not.optval(mask, .true.) .or. size(x) < 2) then
res = ieee_value(1._dp, ieee_quiet_nan)
return
end if
res = 1
end function ${RName}$
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
#:set RName = rname("corr_mask",1, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:)
integer, intent(in) :: dim
logical, intent(in) :: mask(:)
real(${k1}$) :: res
if (count(mask) < 2) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if
res = 1
end function ${RName}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:set RName = rname("corr_mask",1, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:)
integer, intent(in) :: dim
logical, intent(in) :: mask(:)
real(dp) :: res
if (count(mask) < 2) then
res = ieee_value(1._dp, ieee_quiet_nan)
return
end if
res = 1
end function ${RName}$
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
#:set RName = rname("corr",2, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:, :)
integer, intent(in) :: dim
logical, intent(in), optional :: mask
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
, merge(size(x, 1), size(x, 2), mask = 1<dim))
integer :: i, j
${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
${t1}$ :: center(size(x, 1),size(x, 2))
if (.not.optval(mask, .true.) .or. size(x) < 2) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if
mean_ = mean(x, dim)
select case(dim)
case(1)
do i = 1, size(x, 1)
center(i, :) = x(i, :) - mean_
end do
#:if t1[0] == 'r'
res = matmul( transpose(center), center)
#:else
res = matmul( transpose(conjg(center)), center)
#:endif
case(2)
do i = 1, size(x, 2)
center(:, i) = x(:, i) - mean_
end do
#:if t1[0] == 'r'
res = matmul( center, transpose(center))
#:else
res = matmul( center, transpose(conjg(center)))
#:endif
case default
call error_stop("ERROR (corr): wrong dimension")
end select
mean_ = 1 / sqrt(diag(res))
do i = 1, size(res, 1)
do j = 1, size(res, 2)
res(j, i) = res(j, i) * mean_(i) * mean_(j)
end do
end do
end function ${RName}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:set RName = rname("corr",2, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:, :)
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
, merge(size(x, 1), size(x, 2), mask = 1<dim))
integer :: i, j
real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
real(dp) :: center(size(x, 1),size(x, 2))
if (.not.optval(mask, .true.) .or. size(x) < 2) then
res = ieee_value(1._dp, ieee_quiet_nan)
return
end if
mean_ = mean(x, dim)
select case(dim)
case(1)
do i = 1, size(x, 1)
center(i, :) = real(x(i, :), dp) - mean_
end do
res = matmul( transpose(center), center)
case(2)
do i = 1, size(x, 2)
center(:, i) = real(x(:, i), dp) - mean_
end do
res = matmul( center, transpose(center))
case default
call error_stop("ERROR (corr): wrong dimension")
end select
mean_ = 1 / sqrt(diag(res))
do i = 1, size(res, 1)
do j = 1, size(res, 2)
res(j, i) = res(j, i) * mean_(i) * mean_(j)
end do
end do
end function ${RName}$
#:endfor
#:for k1, t1 in RC_KINDS_TYPES
#:set RName = rname("corr_mask",2, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:, :)
integer, intent(in) :: dim
logical, intent(in) :: mask(:,:)
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
, merge(size(x, 1), size(x, 2), mask = 1<dim))
integer :: i, j
${t1}$ :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
${t1}$ :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
select case(dim)
case(1)
do i = 1, size(res, 2)
do j = 1, size(res, 1)
mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0,0,kind=${k1}$),&
#:endif
mask_)
centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0,0,kind=${k1}$),&
#:endif
mask_)
res(j, i) = dot_product( centerj_, centeri_)&
/sqrt(dot_product( centeri_, centeri_)*&
dot_product( centerj_, centerj_))
end do
end do
case(2)
do i = 1, size(res, 2)
do j = 1, size(res, 1)
mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0,0,kind=${k1}$),&
#:endif
mask_)
centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0,0,kind=${k1}$),&
#:endif
mask_)
res(j, i) = dot_product( centeri_, centerj_)&
/sqrt(dot_product( centeri_, centeri_)*&
dot_product( centerj_, centerj_))
end do
end do
case default
call error_stop("ERROR (corr): wrong dimension")
end select
end function ${RName}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
#:set RName = rname("corr_mask",2, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x(:, :)
integer, intent(in) :: dim
logical, intent(in) :: mask(:,:)
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
, merge(size(x, 1), size(x, 2), mask = 1<dim))
integer :: i, j
real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
select case(dim)
case(1)
do i = 1, size(res, 2)
do j = 1, size(res, 1)
mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
0._dp, mask_)
centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
0._dp, mask_)
res(j, i) = dot_product( centerj_, centeri_)&
/sqrt(dot_product( centeri_, centeri_)*&
dot_product( centerj_, centerj_))
end do
end do
case(2)
do i = 1, size(res, 2)
do j = 1, size(res, 1)
mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
0._dp, mask_)
centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
0._dp, mask_)
res(j, i) = dot_product( centeri_, centerj_)&
/sqrt(dot_product( centeri_, centeri_)*&
dot_product( centerj_, centerj_))
end do
end do
case default
call error_stop("ERROR (corr): wrong dimension")
end select
end function ${RName}$
#:endfor
end submodule
