test_median.fypp

#:include "common.fypp"
#:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES

#:set NRANK = 3

module test_stats_median
    use testdrive, only : new_unittest, unittest_type, error_type, check
    use stdlib_stats, only: median
    use stdlib_kinds, only : int8, int16, int32, int64, sp, dp, xdp, qp
    use, intrinsic :: ieee_arithmetic, only : ieee_is_nan
    implicit none
    private

    public :: collect_stats_median

    real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
    real(dp), parameter :: dptol = 2000 * epsilon(1._dp)
#:if WITH_XDP
    real(xdp), parameter :: xdptol = 2000 * epsilon(1._xdp)
#:endif
#:if WITH_QP
    real(qp), parameter :: qptol = 2000 * epsilon(1._qp)
#:endif

    #:for k1,t1 in IR_KINDS_TYPES
    ${t1}$ , parameter :: d1_${k1}$(12) = [${t1}$ :: 10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20]
    ${t1}$ :: d2_${k1}$(3, 4) = reshape(d1_${k1}$, [3, 4])
    ${t1}$ :: d3_${k1}$(2, 3, 2) = reshape(d1_${k1}$, [2, 3, 2])
    ${t1}$ , parameter :: d1odd_${k1}$(13) = [${t1}$ :: d1_${k1}$, 20]
    ${t1}$ :: d2odd_${k1}$(3, 5) = reshape(d1odd_${k1}$, [3, 5], [${t1}$ :: 0])
    ${t1}$ :: d3odd_${k1}$(1, 3, 5) = reshape(d1odd_${k1}$, [1, 3, 5], [${t1}$ :: 0])
    #:endfor


contains

    !> Collect all exported unit tests
    subroutine collect_stats_median(testsuite)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            new_unittest("test_stats_median_size_int8", test_stats_median_size_int8) &
            #:for k1,t1 in IR_KINDS_TYPES
            , new_unittest("test_stats_median_size_${k1}$", test_stats_median_size_${k1}$) &
            , new_unittest("test_stats_median_odd_size_${k1}$", test_stats_median_odd_size_${k1}$) &
            , new_unittest("test_stats_median_all_${k1}$", test_stats_median_all_${k1}$) &
            , new_unittest("test_stats_median_all_odd_${k1}$", test_stats_median_all_odd_${k1}$) &
            , new_unittest("test_stats_median_all_optmask_${k1}$", test_stats_median_all_optmask_${k1}$) &
            , new_unittest("test_stats_median_${k1}$", test_stats_median_${k1}$) &
            , new_unittest("test_stats_median_odd_${k1}$", test_stats_median_odd_${k1}$) &
            , new_unittest("test_stats_median_optmask_${k1}$", test_stats_median_optmask_${k1}$) &
            , new_unittest("test_stats_median_mask_all_${k1}$", test_stats_median_mask_all_${k1}$) &
            , new_unittest("test_stats_median_mask_${k1}$", test_stats_median_mask_${k1}$) &
            #:endfor
            ]
    end subroutine collect_stats_median

    #:for k1,t1 in INT_KINDS_TYPES
    subroutine test_stats_median_size_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))
        !check just to be sure that the setup of d0 is correct
        call check(error, size(d0), 0, 'size(d0): should be of size 0')

        #:for rank in range(1, NRANK + 1)
        call check(error, mod(size(d${rank}$_${k1}$), 2), 0&
                    , 'mod(size(d${rank}$_${k1}$), 2): should be an even number'&
                    )
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_odd_size_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        #:for rank in range(1, NRANK + 1)
        call check(error, mod(size(d${rank}$odd_${k1}$), 2), 1&
                    , 'mod(size(d${rank}$)_${k1}$, 2): should be an odd number'&
                    )
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_all_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))
        call check(error, ieee_is_nan(median(d0)), 'median(d0): should be NaN' )
        if (allocated(error)) return

        #:for rank in range(1, NRANK + 1)
        call check(error, median(d${rank}$_${k1}$),  1.5_dp&
                    , 'median(d${rank}$_${k1}$): uncorrect answer'&
                    , thr = dptol)
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_all_odd_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, median(d1odd_${k1}$),  2._dp&
                    , 'median(d1odd_${k1}$): uncorrect answer'&
                    , thr = dptol)
        if (allocated(error)) return

        call check(error, median(d2odd_${k1}$),  1._dp&
                    , 'median(d2odd_${k1}$): uncorrect answer'&
                    , thr = dptol)

    end subroutine

    subroutine test_stats_median_all_optmask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0_${k1}$(:)

        allocate(d0_${k1}$(0))

        #:for rank in range(0, NRANK + 1)
        call check(error, ieee_is_nan(median(d${rank}$_${k1}$, .false.))&
                    , 'median(d${rank}$_${k1}$, .false.): uncorrect answer')
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error, ieee_is_nan(median(d0, 1)), 'median(d0, 1): should return NaN' )

        call check(error&
                    , abs(median(d1_${k1}$, 1) - 1.5_dp) < dptol&
                    , 'median(d1_${k1}$, 1): uncorrect answer'&
                    )
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2_${k1}$, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol&
                    , 'median(d2_${k1}$, 1): uncorrect answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2_${k1}$, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol&
                    ,'median(d2_${k1}$, 2): uncorrect answer')
        if (allocated(error)) return

    end subroutine

    subroutine test_stats_median_odd_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error&
                    , abs(median(d1odd_${k1}$, 1) - 2._dp) < dptol&
                    , 'median(d1odd_${k1}$, 1): wrong answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2odd_${k1}$, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol&
                    , 'median(d2odd_${k1}$, 1): wrong answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2odd_${k1}$, 2) - [7._dp, 1._dp, 0._dp])) < dptol&
                    , 'median(d2odd_${k1}$, 2): wrong answer')
        if (allocated(error)) return

    end subroutine

    subroutine test_stats_median_optmask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error, ieee_is_nan(median(d0, 1, .false.))&
                    , 'median(d0, 1, .false.): uncorrect answer'&
                    )
        if (allocated(error)) return

        call check(error, ieee_is_nan(median(d1_${k1}$, 1, .false.))&
                    , 'median(d1_${k1}$, 1, .false.): uncorrect answer'&
                    )
        if (allocated(error)) return

        #:for rank in range(2, NRANK+1)
        #:for dim in range(1, rank+1)
        call check(error, any(ieee_is_nan(median(d${rank}$_${k1}$, ${dim}$, .false.)))&
                    , 'median(d${rank}$_${k1}$, ${dim}$, .false.): uncorrect answer')
        if (allocated(error)) return
        #:endfor
        #:endfor
    end subroutine

    subroutine test_stats_median_mask_all_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error,  ieee_is_nan(median(d0, d0 > 0))&
                    , 'median(d0, d0 > 0): should be NaN' )
        if (allocated(error)) return

        #:for rank in range(1, NRANK+1)
        call check(error&
                    , ieee_is_nan(median(d${rank}$_${k1}$, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$)))&
                    , 'median(d${rank}$_${k1}$, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$))' )
        if (allocated(error)) return
        #:endfor

        #:for rank in range(1, NRANK+1)
        call check(error&
                    , (median(d${rank}$_${k1}$, d${rank}$_${k1}$ > 0) - 7._dp) < dptol&
                    , 'median(d${rank}$_${k1}$, d${rank}$_${k1}$> 0)' ) 
        if (allocated(error)) return
        #:endfor

    end subroutine

    subroutine test_stats_median_mask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error&
                   , ieee_is_nan(median(d0, 1, d0 > 0))&
                   , 'median(d0, 1, d0 > 0): uncorrect answer' )
        if (allocated(error)) return

        call check(error&
                    , ieee_is_nan(median(d1_${k1}$, 1, d1_${k1}$ > huge(d1_${k1}$)))&
                    , 'median(d1_${k1}$, 1_${k1}$, d1_${k1}$ > huge(d1_${k1}$)): answer should be IEEE NaN' )
        if (allocated(error)) return

        #:for rank in range(2, NRANK+1)
        call check(error&
                    , any(ieee_is_nan(median(d${rank}$_${k1}$, 1, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$))))&
                    , 'median(d${rank}$_${k1}$, 1_${k1}$, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$)): answer should be IEEE NaN' )
        if (allocated(error)) return
        #:endfor

        call check(error&
                    , (median(d1_${k1}$, 1, d1_${k1}$ > 0) - 7._dp) < dptol&
                    , 'median(d1_${k1}$, 1, d1_${k1}$ >0): uncorrect answer') 
        if (allocated(error)) return

        call check(error&
                    , sum(abs( (median(d2_${k1}$, 1, d2_${k1}$ > 0) - [ 6._dp, 6._dp, 8._dp, 10.5_dp] )  )) &
                    < dptol&
                    , 'median(d2_${k1}$, 1, d2_${k1}$ > 0): uncorrect answer') 
        if (allocated(error)) return

        call check(error&
                    , sum(abs((median(d2_${k1}$, 2, d2_${k1}$ > 0) - [ 8.5_dp, 2._dp, 14.5_dp] )))&
                    < dptol&
                    , 'median(d2_${k1}$, 2, d2_${k1}$ > 0)') 
        if (allocated(error)) return

        call check(error&
                    , any(ieee_is_nan(median(d3_${k1}$, 1, d3_${k1}$ > 0)))&
                    , 'median(d3_${k1}$, 1, d3_${k1}$ > 0): should contain at least 1 IEEE NaN')

    end subroutine
    #:endfor

    #:for k1,t1 in REAL_KINDS_TYPES
    subroutine test_stats_median_size_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))
        !check just to be sure that the setup of d0 is correct
        call check(error, size(d0), 0, 'size(d0): should be of size 0')

        #:for rank in range(1, NRANK + 1)
        call check(error, mod(size(d${rank}$_${k1}$), 2), 0&
                    , 'mod(size(d${rank}$_${k1}$), 2): should be an even number'&
                    )
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_odd_size_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        #:for rank in range(1, NRANK + 1)
        call check(error, mod(size(d${rank}$odd_${k1}$), 2), 1&
                    , 'mod(size(d${rank}$_${k1}$), 2): should be an odd number'&
                    )
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_all_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))
        call check(error,  ieee_is_nan(median(d0)), 'median(d0): should be NaN' )
        if (allocated(error)) return

        #:for rank in range(1, NRANK + 1)
        call check(error, median(d${rank}$_${k1}$),  1.5_${k1}$&
                    , 'median(d${rank}$_${k1}$): uncorrect answer'&
                    , thr = ${k1}$tol)
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_all_odd_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error, median(d1odd_${k1}$),  2._${k1}$&
                    , 'median(d1odd_${k1}$): uncorrect answer'&
                    , thr = ${k1}$tol)
        if (allocated(error)) return

        call check(error, median(d2odd_${k1}$),  1._${k1}$&
                    , 'median(d2odd_${k1}$): uncorrect answer'&
                    , thr = ${k1}$tol)
        if (allocated(error)) return

        call check(error, median(d2odd_${k1}$),  1._${k1}$&
                    , 'median(d2odd_${k1}$): uncorrect answer'&
                    , thr = ${k1}$tol)
        if (allocated(error)) return

    end subroutine

    subroutine test_stats_median_all_optmask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0_${k1}$(:)

        allocate(d0_${k1}$(0))

        #:for rank in range(0, NRANK + 1)
        call check(error, ieee_is_nan(median(d${rank}$_${k1}$, .false.))&
                    , 'median(d${rank}$_${k1}$, .false.): uncorrect answer')
        if (allocated(error)) return
        #:endfor
    end subroutine

    subroutine test_stats_median_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error, ieee_is_nan(median(d0, 1)), 'median(d0, 1): should return NaN' )

        call check(error&
                    , abs(median(d1_${k1}$, 1) - 1.5_${k1}$) < ${k1}$tol&
                    , 'median(d1_${k1}$, 1): uncorrect answer'&
                    )
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2_${k1}$, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$])) < ${k1}$tol&
                    , 'median(d2_${k1}$, 1): uncorrect answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2_${k1}$, 2) - [3.5_${k1}$, 1.5_${k1}$, 3._${k1}$])) < ${k1}$tol&
                    ,'median(d2_${k1}$, 2): uncorrect answer')
        if (allocated(error)) return

    end subroutine

    subroutine test_stats_median_odd_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        call check(error&
                    , abs(median(d1odd_${k1}$, 1) - 2._${k1}$) < ${k1}$tol&
                    , 'median(d1odd_${k1}$, 1): wrong answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2odd_${k1}$, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol&
                    , 'median(d2odd_${k1}$, 1): wrong answer')
        if (allocated(error)) return

        call check(error&
                    , sum(abs(median(d2odd_${k1}$, 2) - [7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol&
                    , 'median(d2odd_${k1}$, 2): wrong answer')
        if (allocated(error)) return

    end subroutine

    subroutine test_stats_median_optmask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error, ieee_is_nan(median(d0, 1, .false.))&
                    , 'median(d0, 1, .false.): uncorrect answer'&
                    )
        if (allocated(error)) return

        call check(error, ieee_is_nan(median(d1_${k1}$, 1, .false.))&
                    , 'median(d1_${k1}$, 1, .false.): uncorrect answer'&
                    )
        if (allocated(error)) return

        #:for rank in range(2, NRANK+1)
        #:for dim in range(1, rank+1)
        call check(error, any(ieee_is_nan(median(d${rank}$_${k1}$, ${dim}$, .false.)))&
                    , 'median(d${rank}$_${k1}$, ${dim}$, .false.): uncorrect answer')
        if (allocated(error)) return
        #:endfor
        #:endfor
    end subroutine

    subroutine test_stats_median_mask_all_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error, ieee_is_nan(median(d0, d0 > 0))&
                    , 'median(d0, d0 > 0): should be NaN' )
        if (allocated(error)) return

        #:for rank in range(1, NRANK+1)
        call check(error&
                    , ieee_is_nan(median(d${rank}$_${k1}$, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$)))&
                    , 'median(d${rank}$_${k1}$, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$))' )
        if (allocated(error)) return
        #:endfor

        #:for rank in range(1, NRANK+1)
        call check(error&
                    , (median(d${rank}$_${k1}$, d${rank}$_${k1}$ > 0) - 7._${k1}$) < ${k1}$tol&
                    , 'median(d${rank}$_${k1}$, d${rank}$_${k1}$ > 0)' ) 
        if (allocated(error)) return
        #:endfor

    end subroutine

    subroutine test_stats_median_mask_${k1}$(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        ${t1}$, allocatable :: d0(:)

        allocate(d0(0))

        call check(error&
                   , ieee_is_nan(median(d0, 1, d0 > 0))&
                   , 'median(d0, 1, d0 > 0): uncorrect answer' )
        if (allocated(error)) return

        call check(error&
                    , ieee_is_nan(median(d1_${k1}$, 1, d1_${k1}$ > huge(d1_${k1}$)))&
                    , 'median(d1_${k1}$, 1, d1_${k1}$ > huge(d1_${k1}$)): answer should be IEEE NaN' )
        if (allocated(error)) return

        #:for rank in range(2, NRANK+1)
        call check(error&
                    , any(ieee_is_nan(median(d${rank}$_${k1}$, 1, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$))))&
                    , 'median(d${rank}$_${k1}$, 1, d${rank}$_${k1}$ > huge(d${rank}$_${k1}$)): answer should be IEEE NaN' )
        if (allocated(error)) return
        #:endfor

        call check(error&
                    , (median(d1_${k1}$, 1, d1_${k1}$ > 0) - 7._${k1}$) < ${k1}$tol&
                    , 'median(d1_${k1}$, 1, d1_${k1}$ >0): uncorrect answer') 
        if (allocated(error)) return

        call check(error&
                    , sum(abs( (median(d2_${k1}$, 1, d2_${k1}$ > 0) - [ 6._${k1}$, 6._${k1}$, 8._${k1}$, 10.5_${k1}$] )  )) &
                    < ${k1}$tol&
                    , 'median(d2_${k1}$, 1, d2_${k1}$ > 0): uncorrect answer') 
        if (allocated(error)) return

        call check(error&
                    , sum(abs((median(d2_${k1}$, 2, d2_${k1}$ > 0) - [ 8.5_${k1}$, 2._${k1}$, 14.5_${k1}$] )))&
                    < ${k1}$tol&
                    , 'median(d2_${k1}$, 2, d2_${k1}$ > 0)') 
        if (allocated(error)) return

        call check(error&
                    , any(ieee_is_nan(median(d3_${k1}$, 1, d3_${k1}$ > 0)))&
                    , 'median(d3_${k1}$, 1, d3_${k1}$ > 0): should contain at least 1 IEEE NaN')

    end subroutine

    #:endfor

end module test_stats_median

program tester
    use, intrinsic :: iso_fortran_env, only : error_unit
    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
    use test_stats_median, only : collect_stats_median
    implicit none
    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("stats_median", collect_stats_median) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if
end program