stdlib_stats_cov.fypp Source File


This file depends on

sourcefile~~stdlib_stats_cov.fypp~~EfferentGraph sourcefile~stdlib_stats_cov.fypp stdlib_stats_cov.fypp sourcefile~stdlib_stats.fypp stdlib_stats.fypp sourcefile~stdlib_stats_cov.fypp->sourcefile~stdlib_stats.fypp sourcefile~stdlib_optval.fypp stdlib_optval.fypp sourcefile~stdlib_stats_cov.fypp->sourcefile~stdlib_optval.fypp sourcefile~stdlib_error.f90 stdlib_error.f90 sourcefile~stdlib_stats_cov.fypp->sourcefile~stdlib_error.f90 sourcefile~stdlib_kinds.f90 stdlib_kinds.f90 sourcefile~stdlib_stats.fypp->sourcefile~stdlib_kinds.f90 sourcefile~stdlib_optval.fypp->sourcefile~stdlib_kinds.f90 sourcefile~stdlib_error.f90->sourcefile~stdlib_optval.fypp

Contents

Source Code


Source Code

#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_stats) stdlib_stats_cov

  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
  use stdlib_error, only: error_stop
  use stdlib_optval, only: optval
  implicit none

contains

  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("cov",1, t1, k1)
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      logical, intent(in), optional :: corrected
      real(${k1}$) :: res

      if (.not.optval(mask, .true.)) then
        res = ieee_value(1._${k1}$, ieee_quiet_nan)
        return
      end if

      res = var(x, dim, corrected = optval(corrected, .true.))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in INT_KINDS_TYPES
    #:set RName = rname("cov",1, t1, k1, 'dp')
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      logical, intent(in), optional :: corrected
      real(dp) :: res

      if (.not.optval(mask, .true.)) then
        res = ieee_value(1._dp, ieee_quiet_nan)
        return
      end if

      res = var(x, dim, corrected = optval(corrected, .true.))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("cov_mask",1, t1, k1)
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      logical, intent(in), optional :: corrected
      real(${k1}$) :: res

      res = var(x, dim, mask, corrected = optval(corrected, .true.))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in INT_KINDS_TYPES
    #:set RName = rname("cov_mask",1, t1, k1, 'dp')
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:)
      logical, intent(in), optional :: corrected
      real(dp) :: res

      res = var(x, dim, mask, corrected = optval(corrected, .true.))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("cov",2, t1, k1)
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      logical, intent(in), optional :: corrected
      ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
                          , merge(size(x, 1), size(x, 2), mask = 1<dim))

      integer :: i
      ${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.)) 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 (cov): wrong dimension")
      end select
      res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in INT_KINDS_TYPES
    #:set RName = rname("cov",2, t1, k1, 'dp')
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in), optional :: mask
      logical, intent(in), optional :: corrected
      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
      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.)) 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 (cov): wrong dimension")
      end select
      res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))

    end function ${RName}$
  #:endfor


  #:for k1, t1 in RC_KINDS_TYPES
    #:set RName = rname("cov_mask",2, t1, k1)
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      logical, intent(in), optional :: corrected
      ${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, n
      ${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_)

              n = count(mask_)
              res(j, i) = dot_product( centerj_, centeri_)&
                           / (n - merge(1, 0,&
                            optval(corrected, .true.) .and. n > 0))
            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_)

              n = count(mask_)
              res(j, i) = dot_product( centeri_, centerj_)&
                           / (n - merge(1, 0,&
                            optval(corrected, .true.) .and. n > 0))
            end do
          end do
        case default
          call error_stop("ERROR (cov): wrong dimension")
      end select

    end function ${RName}$
  #:endfor


  #:for k1, t1 in INT_KINDS_TYPES
    #:set RName = rname("cov_mask",2, t1, k1, 'dp')
    module function ${RName}$(x, dim, mask, corrected) result(res)
      ${t1}$, intent(in) :: x(:, :)
      integer, intent(in) :: dim
      logical, intent(in) :: mask(:,:)
      logical, intent(in), optional :: corrected
      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, n
      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_)

              n = count(mask_)
              res(j, i) = dot_product( centerj_, centeri_)&
                           / (n - merge(1, 0,&
                            optval(corrected, .true.) .and. n > 0))
            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_)

              n = count(mask_)
              res(j, i) = dot_product( centeri_, centerj_)&
                           / (n - merge(1, 0,&
                            optval(corrected, .true.) .and. n > 0))
            end do
          end do
        case default
          call error_stop("ERROR (cov): wrong dimension")
      end select

    end function ${RName}$
  #:endfor


end submodule