From 8dddff48818a962dc6ed124bb7cb2015442e5f7b Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Thu, 3 Jun 2021 23:29:54 +0200 Subject: [PATCH 01/30] progress --- src/CMakeLists.txt | 1 + src/Makefile.manual | 6 + src/common.fypp | 30 +++++ src/stdlib_stats.fypp | 58 +++++++++- src/stdlib_stats_median.fypp | 188 ++++++++++++++++++++++++++++++++ src/tests/stats/CMakeLists.txt | 1 + src/tests/stats/test_median.f90 | 113 +++++++++++++++++++ 7 files changed, 396 insertions(+), 1 deletion(-) create mode 100644 src/stdlib_stats_median.fypp create mode 100644 src/tests/stats/test_median.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fa7f6e639..046d1f7f8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,6 +18,7 @@ set(fppFiles stdlib_stats_corr.fypp stdlib_stats_cov.fypp stdlib_stats_mean.fypp + stdlib_stats_median.fypp stdlib_stats_moment.fypp stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index e55bb80fd..d64682508 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,6 +18,7 @@ SRCFYPP =\ stdlib_stats_corr.fypp \ stdlib_stats_cov.fypp \ stdlib_stats_mean.fypp \ + stdlib_stats_median.fypp \ stdlib_stats_moment.fypp \ stdlib_stats_moment_all.fypp \ stdlib_stats_moment_mask.fypp \ @@ -108,6 +109,11 @@ stdlib_stats_mean.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_median.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_sorting.o \ + stdlib_stats.o stdlib_stats_moment.o: \ stdlib_optval.o \ stdlib_kinds.o \ diff --git a/src/common.fypp b/src/common.fypp index 42be93691..cd4864927 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -147,4 +147,34 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endcall #:enddef + +#:def select_subvector(rank, idim) + #:assert rank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, idim) + ${"j" + "_" * (i) }$ + #:endfor + : + #:for i in range(idim + 1, rank + 1) + ${"j" + "_" * (i) }$ + #:endfor + #:endcall +#:enddef + +#:def reduce_subvector(rank, idim) + #:assert rank > 0 + #:if rank > 1 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, idim) + ${"j" + "_" * (i) }$ + #:endfor + #:for i in range(idim + 1, rank + 1) + ${"j" + "_" * (i) }$ + #:endfor + #:endcall + #:endif +#:enddef + + + #:endmute diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 4aa9edaa1..99f939ba9 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -11,7 +11,7 @@ module stdlib_stats implicit none private ! Public API - public :: corr, cov, mean, moment, var + public :: corr, cov, mean, median, moment, var interface corr @@ -304,6 +304,62 @@ module stdlib_stats end interface mean + interface median + !! version: experimental + !! + !! Median of array elements + !! ([Specification](../page/specs/stdlib_stats.html#median)) + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median_all",rank, t1, k1) + module function ${RName}$ (x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + ${t1}$ :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median_all",rank, t1, k1,'dp') + module function ${RName}$ (x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median",rank, t1, k1,'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + + + + end interface + interface var !! version: experimental diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp new file mode 100644 index 000000000..97fe14da0 --- /dev/null +++ b/src/stdlib_stats_median.fypp @@ -0,0 +1,188 @@ +#:include "common.fypp" +#:set RANKS = range(1, MAXRANK + 1) +submodule (stdlib_stats) stdlib_stats_median + + use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use stdlib_error, only: error_stop + use stdlib_optval, only: optval + use stdlib_sorting, only: sort + implicit none + +contains + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median_all",rank, t1, k1) + module function ${RName}$ (x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + ${t1}$ :: res + + integer :: c, n + ${t1}$, allocatable :: x_tmp(:) + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + x_tmp = reshape(x, [size(x)]) + + n = size(x_tmp) + c = floor( (n + 1) / 2._${k1}$ ) + + call sort(x_tmp) + + if (mod(n, 2) == 0) then + res = sum(x_tmp(c:c+1)) / 2 + else + res = x_tmp(c) + endif + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median_all",rank, t1, k1,'dp') + module function ${RName}$ (x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + + integer :: c, n + ${t1}$, allocatable :: x_tmp(:) + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + x_tmp = reshape(x, [size(x)]) + + n = size(x_tmp) + c = floor( (n + 1) / 2.) + + call sort(x_tmp) + + if (mod(n, 2) == 0) then + res = sum(real(x_tmp(c:c+1), kind=dp)) / 2._dp + else + res = x_tmp(c) + endif + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: c, n + #:for fj in range(1, rank+1) + integer :: j${"_" * fj}$ + #:endfor + ${t1}$, allocatable :: x_tmp(:) + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + n = size(x, dim) + c = floor( (n + 1) / 2._${k1}$ ) + + allocate(x_tmp(n)) + + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + #:for fj in range(1, fi) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + #:for fj in range(fi+1, rank+1) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + x_tmp(:) = x${select_subvector(rank, fi)}$ + call sort(x_tmp) + + if (mod(n, 2) == 0) then + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${k1}$ + else + res${reduce_subvector(rank, fi)}$ = x_tmp(c) + endif + #:for fj in range(1, rank) + end do + #:endfor + #:endfor + case default + call error_stop("ERROR (var): wrong dimension") + end select + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("median",rank, t1, k1,'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: c, n + #:for fj in range(1, rank+1) + integer :: j${"_" * fj}$ + #:endfor + ${t1}$, allocatable :: x_tmp(:) + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + n = size(x, dim) + c = floor( (n + 1) / 2._dp ) + + allocate(x_tmp(n)) + + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + #:for fj in range(1, fi) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + #:for fj in range(fi+1, rank+1) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + x_tmp(:) = x${select_subvector(rank, fi)}$ + call sort(x_tmp) + + if (mod(n, 2) == 0) then + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._dp + else + res${reduce_subvector(rank, fi)}$ = x_tmp(c) + endif + #:for fj in range(1, rank) + end do + #:endfor + #:endfor + case default + call error_stop("ERROR (var): wrong dimension") + end select + + end function ${RName}$ + #:endfor + #:endfor + + +end submodule diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 38f9bb84b..ab7f5085c 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,6 +1,7 @@ ADDTEST(corr) ADDTEST(cov) ADDTEST(mean) +ADDTEST(median) ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) diff --git a/src/tests/stats/test_median.f90 b/src/tests/stats/test_median.f90 new file mode 100644 index 000000000..9cf67efdb --- /dev/null +++ b/src/tests/stats/test_median.f90 @@ -0,0 +1,113 @@ +program test_median + use stdlib_error, only: check + use stdlib_kinds, only: int32, sp, dp + use stdlib_stats, only: median + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 2000 * epsilon(1._dp) + + call test_median_int32() + call test_median_sp() + +contains +subroutine test_median_int32() + integer(int32), allocatable :: d1(:), d2(:,:), d3(:,:,:) + + !even number + d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] + call check(mod(size(d1), 2) == 0, 'should be an even number') + call check( abs(median(d1) - 1.5_dp) < dptol, 'int32 median(d1), even') + + d2 = reshape(d1, [3, 4]) + call check(mod(size(d2), 2) == 0, 'should be an even number') + call check( abs(median(d2) - 1.5_dp) < dptol, 'int32 median(d2), even') + + d3 = reshape(d1, [2, 3, 2]) + call check(mod(size(d3), 2) == 0, 'should be an even number') + call check( abs(median(d3) - 1.5_dp) < dptol, 'int32 median(d3), even') + + !following dimension + call check( abs(median(d1, 1) - 1.5_dp) < dptol, 'int32 median(d1, 1), even') + + call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol, & + 'int32 median(d2, 1)') + call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol, & + 'int32 median(d2, 2)') + + + + + !odd number + d1 = [d1, 20] + call check(mod(size(d1), 2) == 1, 'should be an odd number') + call check( abs(median(d1) - 2._dp) < dptol, 'int32 median(d1), odd') + + d2 = reshape(d1, [3, 5], pad = [0_int32]) + call check(mod(size(d2), 2) == 1, 'should be an odd number') + call check( abs(median(d2) - 1._dp) < dptol, 'int32 median(d2), odd') + + d3 = reshape(d1, [1, 3, 5], pad = [0_int32]) + call check(mod(size(d3), 2) == 1, 'should be an odd number') + call check( abs(median(d3) - 1._dp) < dptol, 'int32 median(d3), odd') + + !following dimension + call check( abs(median(d1, 1) - 2._dp) < dptol, 'int32 median(d1, 1)') + + call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol, & + 'int32 median(d2, 1), odd') + call check( sum(abs(median(d2, 2) - [7._dp, 1._dp, 0._dp])) < dptol, & + 'int32 median(d2, 2), odd') + +end subroutine + +subroutine test_median_sp() + real(sp), allocatable :: d1(:), d2(:,:), d3(:,:,:) + + !even number + d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] + call check(mod(size(d1), 2) == 0, 'should be an even number') + call check( abs(median(d1) - 1.5_sp) < sptol, 'sp median(d1), even') + + d2 = reshape(d1, [3, 4]) + call check(mod(size(d2), 2) == 0, 'should be an even number') + call check( abs(median(d2) - 1.5_sp) < sptol, 'sp median(d2), even') + + d3 = reshape(d1, [2, 3, 2]) + call check(mod(size(d3), 2) == 0, 'should be an even number') + call check( abs(median(d3) - 1.5_sp) < sptol, 'sp median(d3), even') + + + !following dimension + call check( abs(median(d1, 1) - 1.5_sp) < sptol, 'sp median(d1, 1), even') + + call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp])) < sptol, & + 'sp median(d2, 1)') + call check( sum(abs(median(d2, 2) - [3.5_sp, 1.5_sp, 3._sp])) < sptol, & + 'sp median(d2, 2)') + + + !odd number + d1 = [d1, 20.] + call check(mod(size(d1), 2) == 1, 'should be an odd number') + call check( abs(median(d1) - 2._sp) < sptol, 'sp median(d1), odd') + + d2 = reshape(d1, [3, 5], pad = [0._sp]) + call check(mod(size(d2), 2) == 1, 'should be an odd number') + call check( abs(median(d2) - 1._sp) < sptol, 'sp median(d2), odd') + + d3 = reshape(d1, [1, 3, 5], pad = [0._sp]) + call check(mod(size(d3), 2) == 1, 'should be an odd number') + call check( abs(median(d3) - 1._sp) < sptol, 'sp median(d3), odd') + + !following dimension + call check( abs(median(d1, 1) - 2._sp) < sptol, 'sp median(d1, 1)') + + call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp, 0._sp])) < sptol, & + 'sp median(d2, 1), odd') + call check( sum(abs(median(d2, 2) - [7._sp, 1._sp, 0._sp])) < sptol, & + 'sp median(d2, 2), odd') + + +end subroutine + +end program From dbc16af7213837c1a94954fc89073c9ed548bd95 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 4 Jun 2021 23:09:23 +0200 Subject: [PATCH 02/30] median: add subroutine to compute median of elements of an array --- doc/specs/stdlib_stats.md | 89 +++++++++++++++ src/stdlib_stats.fypp | 48 ++++++++ src/stdlib_stats_median.fypp | 194 ++++++++++++++++++++++++++++++-- src/tests/stats/test_median.f90 | 122 ++++++++++++++++---- 4 files changed, 417 insertions(+), 36 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index de475c60d..944a111d4 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -26,6 +26,10 @@ The Pearson correlation between two rows (or columns), say `x` and `y`, of `arra `result = [[stdlib_stats(module):corr(interface)]](array, dim [, mask])` +### Class + +Generic subroutine + ### Arguments `array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`. @@ -83,6 +87,10 @@ The scaling can be changed with the logical argument `corrected`. If `corrected` `result = [[stdlib_stats(module):cov(interface)]](array, dim [, mask [, corrected]])` +### Class + +Generic subroutine + ### Arguments `array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`. @@ -134,6 +142,10 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a `result = [[stdlib_stats(module):mean(interface)]](array, dim [, mask])` +### Class + +Generic subroutine + ### Arguments `array`: Shall be an array of type `integer`, `real`, or `complex`. @@ -166,6 +178,75 @@ program demo_mean end program demo_mean ``` +## `median` - median of array elements + +### Status + +Experimental + +### Description + +Returns the median of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`. + +After that the elements are sorted in an increasing order, e.g. `array_sorted = +sort(array)`, the median of the elements of `array` are defined as, if +`n = size(array)` is an even number: + +``` +median(array) = array_sorted( floor( (n + 1) / 2.)) +``` + +or as if `n` is an odd number: + +``` +median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) ) +``` + +The array is sorted using the subroutine `sort` provided by the `stdlib_sorting` +module. + +### Syntax + +`result = [[stdlib_stats(module):median(interface)]](array [, mask])` + +`result = [[stdlib_stats(module):median(interface)]](array, dim [, mask])` + +### Class + +Generic subroutine + +### Arguments + +`array`: Shall be an array of type `integer` or `real`. + +`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to `n`, where `n` is the rank of `array`. + +`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`. + +### Return value + +If `array` is of type `real`, the result is of the same type as `array`. +If `array` is of type `integer`, the result is of type `real(dp)`. + +If `dim` is absent, a scalar with the median of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned. + +If `mask` is specified, the result is the median of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`. + +### Example + +```fortran +program demo_median + use stdlib_stats, only: median + implicit none + real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ] + real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3]) + print *, median(x) !returns 3.5 + print *, median(y) !returns 3.5 + print *, median(y, 1) !returns [ 1.5, 3.5, 5.5 ] + print *, median(y, 1,y > 3.) !returns [ NaN, 4.0, 5.5 ] +end program demo_median +``` + ## `moment` - central moments of array elements ### Status @@ -199,6 +280,10 @@ The _k_-th order moment about `center` is defined as : `result = [[stdlib_stats(module):moment(interface)]](array, order, dim [, center [, mask]])` +### Class + +Generic subroutine + ### Arguments `array`: Shall be an array of type `integer`, `real`, or `complex`. @@ -264,6 +349,10 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin `result = [[stdlib_stats(module):var(interface)]](array, dim [, mask [, corrected]])` +### Class + +Generic subroutine + ### Arguments `array`: Shall be an array of type `integer`, `real`, or `complex`. diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 99f939ba9..ebf592309 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -355,6 +355,54 @@ module stdlib_stats #:endfor #:endfor + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask_all',rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask_all',rank, t1, k1,'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + end function ${RName}$ + #:endfor + #:endfor + + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask',rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask',rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 97fe14da0..2950a7ca8 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -21,20 +21,20 @@ contains integer :: c, n ${t1}$, allocatable :: x_tmp(:) - if (.not.optval(mask, .true.)) then + if (.not.optval(mask, .true.) .or. size(x) == 0) then res = ieee_value(1._${k1}$, ieee_quiet_nan) return end if x_tmp = reshape(x, [size(x)]) + call sort(x_tmp) + n = size(x_tmp) c = floor( (n + 1) / 2._${k1}$ ) - call sort(x_tmp) - if (mod(n, 2) == 0) then - res = sum(x_tmp(c:c+1)) / 2 + res = sum(x_tmp(c:c+1)) / 2._${k1}$ else res = x_tmp(c) endif @@ -54,18 +54,18 @@ contains integer :: c, n ${t1}$, allocatable :: x_tmp(:) - if (.not.optval(mask, .true.)) then + if (.not.optval(mask, .true.) .or. size(x) == 0) then res = ieee_value(1._dp, ieee_quiet_nan) return end if x_tmp = reshape(x, [size(x)]) - n = size(x_tmp) - c = floor( (n + 1) / 2.) - call sort(x_tmp) + n = size(x_tmp) + c = floor( (n + 1) / 2._dp) + if (mod(n, 2) == 0) then res = sum(real(x_tmp(c:c+1), kind=dp)) / 2._dp else @@ -91,7 +91,7 @@ contains #:endfor ${t1}$, allocatable :: x_tmp(:) - if (.not.optval(mask, .true.)) then + if (.not.optval(mask, .true.) .or. size(x) == 0) then res = ieee_value(1._${k1}$, ieee_quiet_nan) return end if @@ -123,7 +123,7 @@ contains #:endfor #:endfor case default - call error_stop("ERROR (var): wrong dimension") + call error_stop("ERROR (median): wrong dimension") end select end function ${RName}$ @@ -145,7 +145,7 @@ contains #:endfor ${t1}$, allocatable :: x_tmp(:) - if (.not.optval(mask, .true.)) then + if (.not.optval(mask, .true.) .or. size(x) == 0) then res = ieee_value(1._dp, ieee_quiet_nan) return end if @@ -177,7 +177,177 @@ contains #:endfor #:endfor case default - call error_stop("ERROR (var): wrong dimension") + call error_stop("ERROR (median): wrong dimension") + end select + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask_all',rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res + + integer :: c, n + ${t1}$, allocatable :: x_tmp(:) + + x_tmp = pack(x, mask) + + call sort(x_tmp) + + n = size(x_tmp) + c = floor( (n + 1) / 2._${k1}$ ) + + + if (n == 0) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + else if (mod(n, 2) == 0) then + res = sum(x_tmp(c:c+1)) / 2._${k1}$ + else if (mod(n, 2) == 1) then + res = x_tmp(c) + endif + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask_all',rank, t1, k1,'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + + integer :: c, n + ${t1}$, allocatable :: x_tmp(:) + + x_tmp = pack(x, mask) + + call sort(x_tmp) + + n = size(x_tmp) + c = floor( (n + 1) / 2._dp ) + + if (n == 0) then + res = ieee_value(1._dp, ieee_quiet_nan) + else if (mod(n, 2) == 0) then + res = real(sum(x_tmp(c:c+1)), kind = dp) / 2._dp + else if (mod(n, 2) == 1) then + res = x_tmp(c) + endif + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask',rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: c, n + #:for fj in range(1, rank+1) + integer :: j${"_" * fj}$ + #:endfor + ${t1}$, allocatable :: x_tmp(:) + + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + #:for fj in range(1, fi) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + #:for fj in range(fi+1, rank+1) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + x_tmp = pack(x${select_subvector(rank, fi)}$, & + mask${select_subvector(rank, fi)}$) + call sort(x_tmp) + + n = size(x_tmp) + c = floor( (n + 1) / 2._${k1}$ ) + + if (n == 0) then + res${reduce_subvector(rank, fi)}$ = & + ieee_value(1._${k1}$, ieee_quiet_nan) + else if (mod(n, 2) == 0) then + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${k1}$ + else if (mod(n, 2) == 1) then + res${reduce_subvector(rank, fi)}$ = x_tmp(c) + end if + + deallocate(x_tmp) + #:for fj in range(1, rank) + end do + #:endfor + #:endfor + case default + call error_stop("ERROR (median): wrong dimension") + end select + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname('median_mask',rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: c, n + #:for fj in range(1, rank+1) + integer :: j${"_" * fj}$ + #:endfor + ${t1}$, allocatable :: x_tmp(:) + + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + #:for fj in range(1, fi) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + #:for fj in range(fi+1, rank+1) + do j${"_" * fj}$ = 1, size(x, ${fj}$) + #:endfor + x_tmp = pack(x${select_subvector(rank, fi)}$, & + mask${select_subvector(rank, fi)}$) + call sort(x_tmp) + + n = size(x_tmp) + c = floor( (n + 1) / 2._dp ) + + if (n == 0) then + res${reduce_subvector(rank, fi)}$ = & + ieee_value(1._dp, ieee_quiet_nan) + elseif (mod(n, 2) == 0) then + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._dp + elseif (mod(n, 2) == 1) then + res${reduce_subvector(rank, fi)}$ = x_tmp(c) + endif + + deallocate(x_tmp) + #:for fj in range(1, rank) + end do + #:endfor + #:endfor + case default + call error_stop("ERROR (median): wrong dimension") end select end function ${RName}$ diff --git a/src/tests/stats/test_median.f90 b/src/tests/stats/test_median.f90 index 9cf67efdb..cfa4a4e24 100644 --- a/src/tests/stats/test_median.f90 +++ b/src/tests/stats/test_median.f90 @@ -1,7 +1,8 @@ program test_median use stdlib_error, only: check use stdlib_kinds, only: int32, sp, dp - use stdlib_stats, only: median + use stdlib_stats, only: median, mean + use,intrinsic :: ieee_arithmetic, only : ieee_is_nan implicit none real(sp), parameter :: sptol = 1000 * epsilon(1._sp) real(dp), parameter :: dptol = 2000 * epsilon(1._dp) @@ -11,35 +12,72 @@ program test_median contains subroutine test_median_int32() - integer(int32), allocatable :: d1(:), d2(:,:), d3(:,:,:) + integer(int32), allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) + + allocate(d0(0)) + call check(size(d0) == 0, 'should be of size 0') - !even number d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] call check(mod(size(d1), 2) == 0, 'should be an even number') - call check( abs(median(d1) - 1.5_dp) < dptol, 'int32 median(d1), even') d2 = reshape(d1, [3, 4]) call check(mod(size(d2), 2) == 0, 'should be an even number') - call check( abs(median(d2) - 1.5_dp) < dptol, 'int32 median(d2), even') d3 = reshape(d1, [2, 3, 2]) call check(mod(size(d3), 2) == 0, 'should be an even number') + + !median_all + call check( ieee_is_nan(median(d0)), 'int32 median(d0)' ) + call check( ieee_is_nan(median(d0, .false.)), 'int32 median(d0, .false.)' ) + call check( ieee_is_nan(median(d1, .false.)), 'int32 median(d1, .false.)' ) + call check( ieee_is_nan(median(d2, .false.)), 'int32 median(d2, .false.)' ) + call check( ieee_is_nan(median(d3, .false.)), 'int32 median(d3, .false.)' ) + + call check( abs(median(d1) - 1.5_dp) < dptol, 'int32 median(d1), even') + call check( abs(median(d2) - 1.5_dp) < dptol, 'int32 median(d2), even') call check( abs(median(d3) - 1.5_dp) < dptol, 'int32 median(d3), even') - !following dimension - call check( abs(median(d1, 1) - 1.5_dp) < dptol, 'int32 median(d1, 1), even') + !median + call check( ieee_is_nan(median(d0, 1)), 'int32 median(d0, 1)' ) + call check( ieee_is_nan(median(d0, 1, .false.)), 'int32 median(d0, 1, .false.)' ) + call check( ieee_is_nan(median(d1, 1, .false.)), 'int32 median(d1, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 1, .false.))), 'int32 median(d2, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 2, .false.))), 'int32 median(d2, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 1, .false.))), 'int32 median(d3, 1, .false.)' ) + call check( any(ieee_is_nan(median(d3, 2, .false.))), 'int32 median(d3, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 3, .false.))), 'int32 median(d3, 3, .false.)' ) + call check( abs(median(d1, 1) - 1.5_dp) < dptol, 'int32 median(d1, 1), even') call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol, & 'int32 median(d2, 1)') - call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol, & + call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol, & 'int32 median(d2, 2)') - - + !median_mask_all + call check( ieee_is_nan(median(d0, d0 > 0)), 'int32 median(d0, d0 > 0)' ) + call check( ieee_is_nan(median(d1, d1 > huge(d1))), 'int32 median(d1, d1 > huge(d1))' ) + call check( ieee_is_nan(median(d2, d2 > huge(d2))), 'int32 median(d2, d2 > huge(d2))' ) + call check( ieee_is_nan(median(d3, d3 > huge(d3))), 'int32 median(d3, d3 > huge(d3))' ) + call check( (median(d1, d1 > 0) - 7._dp) < dptol, 'int32 median(d1, d1 > 0)' ) + call check( (median(d2, d2 > 0) - 7._dp) < dptol, 'int32 median(d2, d2 > 0)' ) + call check( (median(d3, d3 > 0) - 7._dp) < dptol, 'int32 median(d3, d3 > 0)' ) + + !median mask + call check( ieee_is_nan(median(d0, 1, d0 > 0)), 'int32 median(d0, 1, d0 > 0)' ) + call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), 'int32 median(d1, 1, d1 > huge(d1))' ) + call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), 'int32 median(d2, 1, d2 > huge(d2))' ) + call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), 'int32 median(d3, 1, d3 > huge(d3))' ) + call check( (median(d1, 1, d1 > 0) - 7._dp) < dptol, 'int32 median(d1, 1, d1 >0)') + call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._dp, 6._dp, 8._dp, 10.5_dp] ) )) & + < dptol, 'int32 median(d2, 1, d2 > 0 )') + call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_dp, 2._dp, 14.5_dp] ) )) & + < dptol, 'int32 median(d2, 2, d2 > 0 )') + call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), 'int32 median(d3, 1, d3 > 0)') !odd number d1 = [d1, 20] call check(mod(size(d1), 2) == 1, 'should be an odd number') + call check( abs(median(d1) - 2._dp) < dptol, 'int32 median(d1), odd') d2 = reshape(d1, [3, 5], pad = [0_int32]) @@ -50,9 +88,7 @@ subroutine test_median_int32() call check(mod(size(d3), 2) == 1, 'should be an odd number') call check( abs(median(d3) - 1._dp) < dptol, 'int32 median(d3), odd') - !following dimension call check( abs(median(d1, 1) - 2._dp) < dptol, 'int32 median(d1, 1)') - call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol, & 'int32 median(d2, 1), odd') call check( sum(abs(median(d2, 2) - [7._dp, 1._dp, 0._dp])) < dptol, & @@ -61,34 +97,72 @@ subroutine test_median_int32() end subroutine subroutine test_median_sp() - real(sp), allocatable :: d1(:), d2(:,:), d3(:,:,:) + real(sp), allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) + + allocate(d0(0)) + call check(size(d0) == 0, 'should be of size 0') - !even number d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] call check(mod(size(d1), 2) == 0, 'should be an even number') - call check( abs(median(d1) - 1.5_sp) < sptol, 'sp median(d1), even') d2 = reshape(d1, [3, 4]) call check(mod(size(d2), 2) == 0, 'should be an even number') - call check( abs(median(d2) - 1.5_sp) < sptol, 'sp median(d2), even') d3 = reshape(d1, [2, 3, 2]) call check(mod(size(d3), 2) == 0, 'should be an even number') + + !median_all + call check( ieee_is_nan(median(d0)), 'sp median(d0)' ) + call check( ieee_is_nan(median(d0, .false.)), 'sp median(d0, .false.)' ) + call check( ieee_is_nan(median(d1, .false.)), 'sp median(d1, .false.)' ) + call check( ieee_is_nan(median(d2, .false.)), 'sp median(d2, .false.)' ) + call check( ieee_is_nan(median(d3, .false.)), 'sp median(d3, .false.)' ) + + call check( abs(median(d1) - 1.5_sp) < sptol, 'sp median(d1), even') + call check( abs(median(d2) - 1.5_sp) < sptol, 'sp median(d2), even') call check( abs(median(d3) - 1.5_sp) < sptol, 'sp median(d3), even') + !median + call check( ieee_is_nan(median(d0, 1)), 'sp median(d0, 1)' ) + call check( ieee_is_nan(median(d0, 1, .false.)), 'sp median(d0, 1, .false.)' ) + call check( ieee_is_nan(median(d1, 1, .false.)), 'sp median(d1, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 1, .false.))), 'sp median(d2, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 2, .false.))), 'sp median(d2, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 1, .false.))), 'sp median(d3, 1, .false.)' ) + call check( any(ieee_is_nan(median(d3, 2, .false.))), 'sp median(d3, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 3, .false.))), 'sp median(d3, 3, .false.)' ) - !following dimension call check( abs(median(d1, 1) - 1.5_sp) < sptol, 'sp median(d1, 1), even') - call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp])) < sptol, & 'sp median(d2, 1)') - call check( sum(abs(median(d2, 2) - [3.5_sp, 1.5_sp, 3._sp])) < sptol, & + call check( sum(abs(median(d2, 2) - [3.5_sp, 1.5_sp, 3._sp])) < sptol, & 'sp median(d2, 2)') + !median_mask_all + call check( ieee_is_nan(median(d0, d0 > 0)), 'sp median(d0, d0 > 0)' ) + call check( ieee_is_nan(median(d1, d1 > huge(d1))), 'sp median(d1, d1 > huge(d1))' ) + call check( ieee_is_nan(median(d2, d2 > huge(d2))), 'sp median(d2, d2 > huge(d2))' ) + call check( ieee_is_nan(median(d3, d3 > huge(d3))), 'sp median(d3, d3 > huge(d3))' ) + call check( (median(d1, d1 > 0) - 7._sp) < sptol, 'sp median(d1, d1 > 0)' ) + call check( (median(d2, d2 > 0) - 7._sp) < sptol, 'sp median(d2, d2 > 0)' ) + call check( (median(d3, d3 > 0) - 7._sp) < sptol, 'sp median(d3, d3 > 0)' ) + + !median mask + call check( ieee_is_nan(median(d0, 1, d0 > 0)), 'sp median(d0, 1, d0 > 0)' ) + call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), 'sp median(d1, 1, d1 > huge(d1))' ) + call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), 'sp median(d2, 1, d2 > huge(d2))' ) + call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), 'sp median(d3, 1, d3 > huge(d3))' ) + call check( (median(d1, 1, d1 > 0) - 7._sp) < sptol, 'sp median(d1, 1, d1 >0)') + call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._sp, 6._sp, 8._sp, 10.5_sp] ) )) & + < sptol, 'sp median(d2, 1, d2 > 0 )') + call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_sp, 2._sp, 14.5_sp] ) )) & + < sptol, 'sp median(d2, 2, d2 > 0 )') + call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), 'sp median(d3, 1, d3 > 0)') !odd number - d1 = [d1, 20.] + d1 = [d1, 20._sp] call check(mod(size(d1), 2) == 1, 'should be an odd number') + call check( abs(median(d1) - 2._sp) < sptol, 'sp median(d1), odd') d2 = reshape(d1, [3, 5], pad = [0._sp]) @@ -99,15 +173,15 @@ subroutine test_median_sp() call check(mod(size(d3), 2) == 1, 'should be an odd number') call check( abs(median(d3) - 1._sp) < sptol, 'sp median(d3), odd') - !following dimension call check( abs(median(d1, 1) - 2._sp) < sptol, 'sp median(d1, 1)') - call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp, 0._sp])) < sptol, & 'sp median(d2, 1), odd') - call check( sum(abs(median(d2, 2) - [7._sp, 1._sp, 0._sp])) < sptol, & + call check( sum(abs(median(d2, 2) - [7._sp, 1._sp, 0._sp])) < sptol, & 'sp median(d2, 2), odd') - end subroutine + + + end program From 5d38d8225c611f2b01fa20712c885fd87b6007f0 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 5 Jun 2021 18:51:37 +0200 Subject: [PATCH 03/30] median: mv to ord_sort + combine REAL and INTEGER procedures --- src/stdlib_stats.fypp | 134 +++++++------------- src/stdlib_stats_median.fypp | 233 +++++------------------------------ 2 files changed, 73 insertions(+), 294 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index ebf592309..7c05ef0db 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -2,6 +2,7 @@ #:set RANKS = range(1, MAXRANK + 1) #:set REDRANKS = range(2, MAXRANK + 1) #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set IR_KINDS_TYPES_OUTPUT = list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) + list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) module stdlib_stats !! Provides support for various statistical methods. This includes currently !! descriptive statistics @@ -309,102 +310,51 @@ module stdlib_stats !! !! Median of array elements !! ([Specification](../page/specs/stdlib_stats.html#median)) - #:for k1, t1 in REAL_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1) - module function ${RName}$ (x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in), optional :: mask - ${t1}$ :: res - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1,'dp') - module function ${RName}$ (x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in), optional :: mask - real(dp) :: res - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1) - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in), optional :: mask - ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1,'dp') - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in), optional :: mask - real(dp) :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1) - module function ${RName}$(x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in) :: mask${ranksuffix(rank)}$ - ${t1}$ :: res - end function ${RName}$ + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT + #:for rank in RANKS + #:set RName = rname("median_all",rank, t1, k1, o1) + module function ${RName}$ (x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(${o1}$) :: res + end function ${RName}$ + #:endfor #:endfor - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1,'dp') - module function ${RName}$(x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in) :: mask${ranksuffix(rank)}$ - real(dp) :: res - end function ${RName}$ + + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT + #:for rank in RANKS + #:set RName = rname("median",rank, t1, k1, o1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor #:endfor - #:endfor - - - - #:for k1, t1 in REAL_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1) - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in) :: mask${ranksuffix(rank)}$ - ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ + + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT + #:for rank in RANKS + #:set RName = rname('median_mask_all',rank, t1, k1, o1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${o1}$) :: res + end function ${RName}$ + #:endfor #:endfor - #:endfor - - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in) :: mask${ranksuffix(rank)}$ - real(dp) :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ + + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT + #:for rank in RANKS + #:set RName = rname('median_mask',rank, t1, k1, o1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor #:endfor - #:endfor - - end interface diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 2950a7ca8..0003904c6 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -1,28 +1,31 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) +#:set IR_KINDS_TYPES_OUTPUT = list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) + list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) + + submodule (stdlib_stats) stdlib_stats_median use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan use stdlib_error, only: error_stop use stdlib_optval, only: optval - use stdlib_sorting, only: sort + use stdlib_sorting, only: sort => ord_sort implicit none contains - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1) + #:set RName = rname("median_all",rank, t1, k1, o1) module function ${RName}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask - ${t1}$ :: res + real(${o1}$) :: res integer :: c, n ${t1}$, allocatable :: x_tmp(:) if (.not.optval(mask, .true.) .or. size(x) == 0) then - res = ieee_value(1._${k1}$, ieee_quiet_nan) + res = ieee_value(1._${o1}$, ieee_quiet_nan) return end if @@ -31,113 +34,26 @@ contains call sort(x_tmp) n = size(x_tmp) - c = floor( (n + 1) / 2._${k1}$ ) + c = floor( (n + 1) / 2._${o1}$ ) if (mod(n, 2) == 0) then - res = sum(x_tmp(c:c+1)) / 2._${k1}$ + res = sum(x_tmp(c:c+1)) / 2._${o1}$ else res = x_tmp(c) - endif - - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1,'dp') - module function ${RName}$ (x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in), optional :: mask - real(dp) :: res - - integer :: c, n - ${t1}$, allocatable :: x_tmp(:) - - if (.not.optval(mask, .true.) .or. size(x) == 0) then - res = ieee_value(1._dp, ieee_quiet_nan) - return - end if - - x_tmp = reshape(x, [size(x)]) - - call sort(x_tmp) - - n = size(x_tmp) - c = floor( (n + 1) / 2._dp) - - if (mod(n, 2) == 0) then - res = sum(real(x_tmp(c:c+1), kind=dp)) / 2._dp - else - res = x_tmp(c) - endif - - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1) - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in), optional :: mask - ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ - - integer :: c, n - #:for fj in range(1, rank+1) - integer :: j${"_" * fj}$ - #:endfor - ${t1}$, allocatable :: x_tmp(:) - - if (.not.optval(mask, .true.) .or. size(x) == 0) then - res = ieee_value(1._${k1}$, ieee_quiet_nan) - return end if - n = size(x, dim) - c = floor( (n + 1) / 2._${k1}$ ) - - allocate(x_tmp(n)) - - select case(dim) - #:for fi in range(1, rank+1) - case(${fi}$) - #:for fj in range(1, fi) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor - #:for fj in range(fi+1, rank+1) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor - x_tmp(:) = x${select_subvector(rank, fi)}$ - call sort(x_tmp) - - if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${k1}$ - else - res${reduce_subvector(rank, fi)}$ = x_tmp(c) - endif - #:for fj in range(1, rank) - end do - #:endfor - #:endfor - case default - call error_stop("ERROR (median): wrong dimension") - end select - end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in INT_KINDS_TYPES + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1,'dp') + #:set RName = rname("median",rank, t1, k1, o1) module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask - real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ integer :: c, n #:for fj in range(1, rank+1) @@ -146,12 +62,12 @@ contains ${t1}$, allocatable :: x_tmp(:) if (.not.optval(mask, .true.) .or. size(x) == 0) then - res = ieee_value(1._dp, ieee_quiet_nan) + res = ieee_value(1._${o1}$, ieee_quiet_nan) return end if n = size(x, dim) - c = floor( (n + 1) / 2._dp ) + c = floor( (n + 1) / 2._${o1}$ ) allocate(x_tmp(n)) @@ -168,10 +84,10 @@ contains call sort(x_tmp) if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._dp + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$ else res${reduce_subvector(rank, fi)}$ = x_tmp(c) - endif + end if #:for fj in range(1, rank) end do #:endfor @@ -185,13 +101,13 @@ contains #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1) + #:set RName = rname('median_mask_all',rank, t1, k1, o1) module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ - ${t1}$ :: res + real(${o1}$) :: res integer :: c, n ${t1}$, allocatable :: x_tmp(:) @@ -201,60 +117,28 @@ contains call sort(x_tmp) n = size(x_tmp) - c = floor( (n + 1) / 2._${k1}$ ) - + c = floor( (n + 1) / 2._${o1}$ ) if (n == 0) then - res = ieee_value(1._${k1}$, ieee_quiet_nan) + res = ieee_value(1._${o1}$, ieee_quiet_nan) else if (mod(n, 2) == 0) then - res = sum(x_tmp(c:c+1)) / 2._${k1}$ + res = sum(x_tmp(c:c+1)) / 2._${o1}$ else if (mod(n, 2) == 1) then res = x_tmp(c) - endif - - end function ${RName}$ - #:endfor - #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1,'dp') - module function ${RName}$(x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in) :: mask${ranksuffix(rank)}$ - real(dp) :: res - - integer :: c, n - ${t1}$, allocatable :: x_tmp(:) - - x_tmp = pack(x, mask) - - call sort(x_tmp) - - n = size(x_tmp) - c = floor( (n + 1) / 2._dp ) - - if (n == 0) then - res = ieee_value(1._dp, ieee_quiet_nan) - else if (mod(n, 2) == 0) then - res = real(sum(x_tmp(c:c+1)), kind = dp) / 2._dp - else if (mod(n, 2) == 1) then - res = x_tmp(c) - endif + end if end function ${RName}$ #:endfor #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1) - module function ${RName}$(x, dim, mask) result(res) + #:set RName = rname('median_mask',rank, t1, k1, o1) + module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ - ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ integer :: c, n #:for fj in range(1, rank+1) @@ -276,13 +160,13 @@ contains call sort(x_tmp) n = size(x_tmp) - c = floor( (n + 1) / 2._${k1}$ ) + c = floor( (n + 1) / 2._${o1}$ ) if (n == 0) then res${reduce_subvector(rank, fi)}$ = & - ieee_value(1._${k1}$, ieee_quiet_nan) + ieee_value(1._${o1}$, ieee_quiet_nan) else if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${k1}$ + res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$ else if (mod(n, 2) == 1) then res${reduce_subvector(rank, fi)}$ = x_tmp(c) end if @@ -300,59 +184,4 @@ contains #:endfor #:endfor - - #:for k1, t1 in INT_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - integer, intent(in) :: dim - logical, intent(in) :: mask${ranksuffix(rank)}$ - real(dp) :: res${reduced_shape('x', rank, 'dim')}$ - - integer :: c, n - #:for fj in range(1, rank+1) - integer :: j${"_" * fj}$ - #:endfor - ${t1}$, allocatable :: x_tmp(:) - - select case(dim) - #:for fi in range(1, rank+1) - case(${fi}$) - #:for fj in range(1, fi) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor - #:for fj in range(fi+1, rank+1) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor - x_tmp = pack(x${select_subvector(rank, fi)}$, & - mask${select_subvector(rank, fi)}$) - call sort(x_tmp) - - n = size(x_tmp) - c = floor( (n + 1) / 2._dp ) - - if (n == 0) then - res${reduce_subvector(rank, fi)}$ = & - ieee_value(1._dp, ieee_quiet_nan) - elseif (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._dp - elseif (mod(n, 2) == 1) then - res${reduce_subvector(rank, fi)}$ = x_tmp(c) - endif - - deallocate(x_tmp) - #:for fj in range(1, rank) - end do - #:endfor - #:endfor - case default - call error_stop("ERROR (median): wrong dimension") - end select - - end function ${RName}$ - #:endfor - #:endfor - - end submodule From 935949eb767296fbe6c4cf6d4f4e9344e8e27401 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 5 Jun 2021 19:56:40 +0200 Subject: [PATCH 04/30] median: mv to sort due to issue #428 --- src/stdlib_stats_median.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 0003904c6..b4c4c81cf 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -8,7 +8,7 @@ submodule (stdlib_stats) stdlib_stats_median use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan use stdlib_error, only: error_stop use stdlib_optval, only: optval - use stdlib_sorting, only: sort => ord_sort + use stdlib_sorting, only: sort => sort implicit none contains From d506f4a23b4ce8be8c7460915bfd71e7849b86c1 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 6 Jun 2021 20:10:52 +0200 Subject: [PATCH 05/30] median: some cleaning --- src/stdlib_stats_median.fypp | 60 ++++++++++++++++++++++++------------ 1 file changed, 41 insertions(+), 19 deletions(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index b4c4c81cf..07896b783 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -21,23 +21,27 @@ contains logical, intent(in), optional :: mask real(${o1}$) :: res - integer :: c, n + integer(kind = int64) :: c, n ${t1}$, allocatable :: x_tmp(:) if (.not.optval(mask, .true.) .or. size(x) == 0) then res = ieee_value(1._${o1}$, ieee_quiet_nan) return end if + + n = size(x, kind=int64) + c = floor( (n + 1) / 2._${o1}$, kind=int64 ) - x_tmp = reshape(x, [size(x)]) + x_tmp = reshape(x, [n]) call sort(x_tmp) - n = size(x_tmp) - c = floor( (n + 1) / 2._${o1}$ ) - - if (mod(n, 2) == 0) then + if (mod(n, 2_int64) == 0) then + #:if t1[0] == 'r' res = sum(x_tmp(c:c+1)) / 2._${o1}$ + #:else + res = sum( real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$ + #:endif else res = x_tmp(c) end if @@ -56,9 +60,11 @@ contains real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ integer :: c, n + #:if rank > 1 #:for fj in range(1, rank+1) integer :: j${"_" * fj}$ #:endfor + #:endif ${t1}$, allocatable :: x_tmp(:) if (.not.optval(mask, .true.) .or. size(x) == 0) then @@ -84,7 +90,12 @@ contains call sort(x_tmp) if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$ + res${reduce_subvector(rank, fi)}$ = & + #:if t1[0] == 'r' + sum(x_tmp(c:c+1)) / 2._${o1}$ + #:else + sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$ + #:endif else res${reduce_subvector(rank, fi)}$ = x_tmp(c) end if @@ -109,21 +120,25 @@ contains logical, intent(in) :: mask${ranksuffix(rank)}$ real(${o1}$) :: res - integer :: c, n - ${t1}$, allocatable :: x_tmp(:) + integer(kind = int64) :: c, n + ${t1}$, allocatable :: x_tmp(:) x_tmp = pack(x, mask) call sort(x_tmp) - n = size(x_tmp) - c = floor( (n + 1) / 2._${o1}$ ) + n = size(x_tmp, kind=int64 ) + c = floor( (n + 1) / 2._${o1}$, kind=int64 ) if (n == 0) then res = ieee_value(1._${o1}$, ieee_quiet_nan) - else if (mod(n, 2) == 0) then + else if (mod(n, 2_int64) == 0) then + #:if t1[0] == 'r' res = sum(x_tmp(c:c+1)) / 2._${o1}$ - else if (mod(n, 2) == 1) then + #:else + res = sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$ + #:endif + else if (mod(n, 2_int64) == 1) then res = x_tmp(c) end if @@ -140,10 +155,12 @@ contains logical, intent(in) :: mask${ranksuffix(rank)}$ real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ - integer :: c, n + integer(kind = int64) :: c, n + #:if rank > 1 #:for fj in range(1, rank+1) integer :: j${"_" * fj}$ #:endfor + #:endif ${t1}$, allocatable :: x_tmp(:) select case(dim) @@ -159,15 +176,20 @@ contains mask${select_subvector(rank, fi)}$) call sort(x_tmp) - n = size(x_tmp) - c = floor( (n + 1) / 2._${o1}$ ) + n = size(x_tmp, kind=int64) + c = floor( (n + 1) / 2._${o1}$, kind=int64 ) if (n == 0) then res${reduce_subvector(rank, fi)}$ = & ieee_value(1._${o1}$, ieee_quiet_nan) - else if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${o1}$ - else if (mod(n, 2) == 1) then + else if (mod(n, 2_int64) == 0) then + res${reduce_subvector(rank, fi)}$ = & + #:if t1[0] == 'r' + sum(x_tmp(c:c+1)) / 2._${o1}$ + #:else + sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$ + #:endif + else if (mod(n, 2_int64) == 1) then res${reduce_subvector(rank, fi)}$ = x_tmp(c) end if From 2bc1a8b849b002116ed766bb80f57054874be3d7 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 6 Jun 2021 20:53:56 +0200 Subject: [PATCH 06/30] median: remove trailing whitespaces --- src/stdlib_stats_median.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 07896b783..567bff01f 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -20,7 +20,7 @@ contains ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res - + integer(kind = int64) :: c, n ${t1}$, allocatable :: x_tmp(:) @@ -28,7 +28,7 @@ contains res = ieee_value(1._${o1}$, ieee_quiet_nan) return end if - + n = size(x, kind=int64) c = floor( (n + 1) / 2._${o1}$, kind=int64 ) From ac1a2a2810ed11cd5d5340eda089346b49b99d1b Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 19:50:46 +0200 Subject: [PATCH 07/30] median: add check on shapes between mask and x --- src/stdlib_stats_median.fypp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 567bff01f..a133b335c 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -123,6 +123,11 @@ contains integer(kind = int64) :: c, n ${t1}$, allocatable :: x_tmp(:) + if (any(shape(x) .ne. shape(mask))) then + call error_stop("ERROR (median): shapes of x and mask are different") + error stop + end if + x_tmp = pack(x, mask) call sort(x_tmp) @@ -163,6 +168,11 @@ contains #:endif ${t1}$, allocatable :: x_tmp(:) + if (any(shape(x) .ne. shape(mask))) then + call error_stop("ERROR (median): shapes of x and mask are different") + error stop + end if + select case(dim) #:for fi in range(1, rank+1) case(${fi}$) From fdfd150daebb0f79b4e851dbeada0d093495c5b6 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 19:55:04 +0200 Subject: [PATCH 08/30] median: add pure statement --- src/stdlib_stats.fypp | 2 +- src/stdlib_stats_median.fypp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 7c05ef0db..5a3ed278b 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -313,7 +313,7 @@ module stdlib_stats #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS #:set RName = rname("median_all",rank, t1, k1, o1) - module function ${RName}$ (x, mask) result(res) + pure module function ${RName}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index a133b335c..6d9b6ecc8 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -16,7 +16,7 @@ contains #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS #:set RName = rname("median_all",rank, t1, k1, o1) - module function ${RName}$ (x, mask) result(res) + pure module function ${RName}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res From b19c5379dae0c0bcc7fad846239816ca9a8e4ba4 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 20:10:47 +0200 Subject: [PATCH 09/30] Update src/stdlib_stats_median.fypp --- src/stdlib_stats_median.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 6d9b6ecc8..58a9e79ed 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -125,7 +125,6 @@ contains if (any(shape(x) .ne. shape(mask))) then call error_stop("ERROR (median): shapes of x and mask are different") - error stop end if x_tmp = pack(x, mask) From 8ed99fef6cf9cd8cb5c8837b44bf03c028795dff Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 20:10:58 +0200 Subject: [PATCH 10/30] Update src/stdlib_stats_median.fypp --- src/stdlib_stats_median.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 58a9e79ed..d7ec66af9 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -169,7 +169,6 @@ contains if (any(shape(x) .ne. shape(mask))) then call error_stop("ERROR (median): shapes of x and mask are different") - error stop end if select case(dim) From 0d46361039f47d25c71cd29fe991bfad265e0bbd Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 21:09:47 +0200 Subject: [PATCH 11/30] median: update test --- src/tests/stats/CMakeLists.txt | 18 +++ src/tests/stats/Makefile.manual | 11 +- src/tests/stats/common.fypp | 37 ++++++ src/tests/stats/test_median.f90 | 187 ----------------------------- src/tests/stats/test_median.fypp | 196 +++++++++++++++++++++++++++++++ 5 files changed, 260 insertions(+), 189 deletions(-) create mode 100644 src/tests/stats/common.fypp delete mode 100644 src/tests/stats/test_median.f90 create mode 100644 src/tests/stats/test_median.fypp diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index ab7f5085c..600e925e5 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,3 +1,21 @@ +### Pre-process: .fpp -> .f90 via Fypp + +# Create a list of the files to be preprocessed +set(fppFiles + test_median.fypp +) + +# Custom preprocessor flags +if(DEFINED CMAKE_MAXIMUM_RANK) + set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}") +elseif(f03rank) + set(fyppFlags) +else() + set(fyppFlags "-DVERSION90") +endif() + +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) + ADDTEST(corr) ADDTEST(cov) ADDTEST(mean) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index a0731e8f4..55a21c5f4 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,5 +1,12 @@ +SRCFYPP =\ + test_median.fypp -PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \ +SRCGEN = $(SRCFYPP:.fypp=.f90) + +$(SRCGEN): %.f90: %.fypp common.fypp + fypp $(FYPPFLAGS) $< $@ + +PROGS_SRC = $(SRCGEN) test_mean.f90 test_moment.f90 test_var.f90 \ test_distribution_PRNG.f90 -include ../Makefile.manual.test.mk \ No newline at end of file +include ../Makefile.manual.test.mk diff --git a/src/tests/stats/common.fypp b/src/tests/stats/common.fypp new file mode 100644 index 000000000..9ebaac6b9 --- /dev/null +++ b/src/tests/stats/common.fypp @@ -0,0 +1,37 @@ +#:mute + +#! Real kinds to be considered during templating +#:set REAL_KINDS = ["sp", "dp", "qp"] + +#! Real types to be considered during templating +#:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS] + +#! Collected (kind, type) tuples for real types +#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES)) + +#! Integer kinds to be considered during templating +#:set INT_KINDS = ["int8", "int16", "int32", "int64"] + +#! Integer types to be considered during templating +#:set INT_TYPES = ["integer({})".format(k) for k in INT_KINDS] + +#! Collected (kind, type) tuples for integer types +#:set INT_KINDS_TYPES = list(zip(INT_KINDS, INT_TYPES)) + + +#! Whether Fortran 90 compatible code should be generated +#:set VERSION90 = defined('VERSION90') + + +#! Ranks to be generated when templates are created +#:if not defined('MAXRANK') + #:if VERSION90 + #:set MAXRANK = 7 + #:else + #:set MAXRANK = 15 + #:endif +#:endif + + + +#:endmute diff --git a/src/tests/stats/test_median.f90 b/src/tests/stats/test_median.f90 deleted file mode 100644 index cfa4a4e24..000000000 --- a/src/tests/stats/test_median.f90 +++ /dev/null @@ -1,187 +0,0 @@ -program test_median - use stdlib_error, only: check - use stdlib_kinds, only: int32, sp, dp - use stdlib_stats, only: median, mean - use,intrinsic :: ieee_arithmetic, only : ieee_is_nan - implicit none - real(sp), parameter :: sptol = 1000 * epsilon(1._sp) - real(dp), parameter :: dptol = 2000 * epsilon(1._dp) - - call test_median_int32() - call test_median_sp() - -contains -subroutine test_median_int32() - integer(int32), allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) - - allocate(d0(0)) - call check(size(d0) == 0, 'should be of size 0') - - d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] - call check(mod(size(d1), 2) == 0, 'should be an even number') - - d2 = reshape(d1, [3, 4]) - call check(mod(size(d2), 2) == 0, 'should be an even number') - - d3 = reshape(d1, [2, 3, 2]) - call check(mod(size(d3), 2) == 0, 'should be an even number') - - !median_all - call check( ieee_is_nan(median(d0)), 'int32 median(d0)' ) - call check( ieee_is_nan(median(d0, .false.)), 'int32 median(d0, .false.)' ) - call check( ieee_is_nan(median(d1, .false.)), 'int32 median(d1, .false.)' ) - call check( ieee_is_nan(median(d2, .false.)), 'int32 median(d2, .false.)' ) - call check( ieee_is_nan(median(d3, .false.)), 'int32 median(d3, .false.)' ) - - call check( abs(median(d1) - 1.5_dp) < dptol, 'int32 median(d1), even') - call check( abs(median(d2) - 1.5_dp) < dptol, 'int32 median(d2), even') - call check( abs(median(d3) - 1.5_dp) < dptol, 'int32 median(d3), even') - - !median - call check( ieee_is_nan(median(d0, 1)), 'int32 median(d0, 1)' ) - call check( ieee_is_nan(median(d0, 1, .false.)), 'int32 median(d0, 1, .false.)' ) - call check( ieee_is_nan(median(d1, 1, .false.)), 'int32 median(d1, 1, .false.)' ) - call check( any(ieee_is_nan(median(d2, 1, .false.))), 'int32 median(d2, 1, .false.)' ) - call check( any(ieee_is_nan(median(d2, 2, .false.))), 'int32 median(d2, 2, .false.)' ) - call check( any(ieee_is_nan(median(d3, 1, .false.))), 'int32 median(d3, 1, .false.)' ) - call check( any(ieee_is_nan(median(d3, 2, .false.))), 'int32 median(d3, 2, .false.)' ) - call check( any(ieee_is_nan(median(d3, 3, .false.))), 'int32 median(d3, 3, .false.)' ) - - call check( abs(median(d1, 1) - 1.5_dp) < dptol, 'int32 median(d1, 1), even') - call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol, & - 'int32 median(d2, 1)') - call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol, & - 'int32 median(d2, 2)') - - !median_mask_all - call check( ieee_is_nan(median(d0, d0 > 0)), 'int32 median(d0, d0 > 0)' ) - call check( ieee_is_nan(median(d1, d1 > huge(d1))), 'int32 median(d1, d1 > huge(d1))' ) - call check( ieee_is_nan(median(d2, d2 > huge(d2))), 'int32 median(d2, d2 > huge(d2))' ) - call check( ieee_is_nan(median(d3, d3 > huge(d3))), 'int32 median(d3, d3 > huge(d3))' ) - call check( (median(d1, d1 > 0) - 7._dp) < dptol, 'int32 median(d1, d1 > 0)' ) - call check( (median(d2, d2 > 0) - 7._dp) < dptol, 'int32 median(d2, d2 > 0)' ) - call check( (median(d3, d3 > 0) - 7._dp) < dptol, 'int32 median(d3, d3 > 0)' ) - - !median mask - call check( ieee_is_nan(median(d0, 1, d0 > 0)), 'int32 median(d0, 1, d0 > 0)' ) - call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), 'int32 median(d1, 1, d1 > huge(d1))' ) - call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), 'int32 median(d2, 1, d2 > huge(d2))' ) - call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), 'int32 median(d3, 1, d3 > huge(d3))' ) - call check( (median(d1, 1, d1 > 0) - 7._dp) < dptol, 'int32 median(d1, 1, d1 >0)') - call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._dp, 6._dp, 8._dp, 10.5_dp] ) )) & - < dptol, 'int32 median(d2, 1, d2 > 0 )') - call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_dp, 2._dp, 14.5_dp] ) )) & - < dptol, 'int32 median(d2, 2, d2 > 0 )') - call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), 'int32 median(d3, 1, d3 > 0)') - - !odd number - d1 = [d1, 20] - call check(mod(size(d1), 2) == 1, 'should be an odd number') - - call check( abs(median(d1) - 2._dp) < dptol, 'int32 median(d1), odd') - - d2 = reshape(d1, [3, 5], pad = [0_int32]) - call check(mod(size(d2), 2) == 1, 'should be an odd number') - call check( abs(median(d2) - 1._dp) < dptol, 'int32 median(d2), odd') - - d3 = reshape(d1, [1, 3, 5], pad = [0_int32]) - call check(mod(size(d3), 2) == 1, 'should be an odd number') - call check( abs(median(d3) - 1._dp) < dptol, 'int32 median(d3), odd') - - call check( abs(median(d1, 1) - 2._dp) < dptol, 'int32 median(d1, 1)') - call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol, & - 'int32 median(d2, 1), odd') - call check( sum(abs(median(d2, 2) - [7._dp, 1._dp, 0._dp])) < dptol, & - 'int32 median(d2, 2), odd') - -end subroutine - -subroutine test_median_sp() - real(sp), allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) - - allocate(d0(0)) - call check(size(d0) == 0, 'should be of size 0') - - d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] - call check(mod(size(d1), 2) == 0, 'should be an even number') - - d2 = reshape(d1, [3, 4]) - call check(mod(size(d2), 2) == 0, 'should be an even number') - - d3 = reshape(d1, [2, 3, 2]) - call check(mod(size(d3), 2) == 0, 'should be an even number') - - !median_all - call check( ieee_is_nan(median(d0)), 'sp median(d0)' ) - call check( ieee_is_nan(median(d0, .false.)), 'sp median(d0, .false.)' ) - call check( ieee_is_nan(median(d1, .false.)), 'sp median(d1, .false.)' ) - call check( ieee_is_nan(median(d2, .false.)), 'sp median(d2, .false.)' ) - call check( ieee_is_nan(median(d3, .false.)), 'sp median(d3, .false.)' ) - - call check( abs(median(d1) - 1.5_sp) < sptol, 'sp median(d1), even') - call check( abs(median(d2) - 1.5_sp) < sptol, 'sp median(d2), even') - call check( abs(median(d3) - 1.5_sp) < sptol, 'sp median(d3), even') - - !median - call check( ieee_is_nan(median(d0, 1)), 'sp median(d0, 1)' ) - call check( ieee_is_nan(median(d0, 1, .false.)), 'sp median(d0, 1, .false.)' ) - call check( ieee_is_nan(median(d1, 1, .false.)), 'sp median(d1, 1, .false.)' ) - call check( any(ieee_is_nan(median(d2, 1, .false.))), 'sp median(d2, 1, .false.)' ) - call check( any(ieee_is_nan(median(d2, 2, .false.))), 'sp median(d2, 2, .false.)' ) - call check( any(ieee_is_nan(median(d3, 1, .false.))), 'sp median(d3, 1, .false.)' ) - call check( any(ieee_is_nan(median(d3, 2, .false.))), 'sp median(d3, 2, .false.)' ) - call check( any(ieee_is_nan(median(d3, 3, .false.))), 'sp median(d3, 3, .false.)' ) - - call check( abs(median(d1, 1) - 1.5_sp) < sptol, 'sp median(d1, 1), even') - call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp])) < sptol, & - 'sp median(d2, 1)') - call check( sum(abs(median(d2, 2) - [3.5_sp, 1.5_sp, 3._sp])) < sptol, & - 'sp median(d2, 2)') - - !median_mask_all - call check( ieee_is_nan(median(d0, d0 > 0)), 'sp median(d0, d0 > 0)' ) - call check( ieee_is_nan(median(d1, d1 > huge(d1))), 'sp median(d1, d1 > huge(d1))' ) - call check( ieee_is_nan(median(d2, d2 > huge(d2))), 'sp median(d2, d2 > huge(d2))' ) - call check( ieee_is_nan(median(d3, d3 > huge(d3))), 'sp median(d3, d3 > huge(d3))' ) - call check( (median(d1, d1 > 0) - 7._sp) < sptol, 'sp median(d1, d1 > 0)' ) - call check( (median(d2, d2 > 0) - 7._sp) < sptol, 'sp median(d2, d2 > 0)' ) - call check( (median(d3, d3 > 0) - 7._sp) < sptol, 'sp median(d3, d3 > 0)' ) - - !median mask - call check( ieee_is_nan(median(d0, 1, d0 > 0)), 'sp median(d0, 1, d0 > 0)' ) - call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), 'sp median(d1, 1, d1 > huge(d1))' ) - call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), 'sp median(d2, 1, d2 > huge(d2))' ) - call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), 'sp median(d3, 1, d3 > huge(d3))' ) - call check( (median(d1, 1, d1 > 0) - 7._sp) < sptol, 'sp median(d1, 1, d1 >0)') - call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._sp, 6._sp, 8._sp, 10.5_sp] ) )) & - < sptol, 'sp median(d2, 1, d2 > 0 )') - call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_sp, 2._sp, 14.5_sp] ) )) & - < sptol, 'sp median(d2, 2, d2 > 0 )') - call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), 'sp median(d3, 1, d3 > 0)') - - !odd number - d1 = [d1, 20._sp] - call check(mod(size(d1), 2) == 1, 'should be an odd number') - - call check( abs(median(d1) - 2._sp) < sptol, 'sp median(d1), odd') - - d2 = reshape(d1, [3, 5], pad = [0._sp]) - call check(mod(size(d2), 2) == 1, 'should be an odd number') - call check( abs(median(d2) - 1._sp) < sptol, 'sp median(d2), odd') - - d3 = reshape(d1, [1, 3, 5], pad = [0._sp]) - call check(mod(size(d3), 2) == 1, 'should be an odd number') - call check( abs(median(d3) - 1._sp) < sptol, 'sp median(d3), odd') - - call check( abs(median(d1, 1) - 2._sp) < sptol, 'sp median(d1, 1)') - call check( sum(abs(median(d2, 1) - [2._sp, -4._sp, 7._sp, 1._sp, 0._sp])) < sptol, & - 'sp median(d2, 1), odd') - call check( sum(abs(median(d2, 2) - [7._sp, 1._sp, 0._sp])) < sptol, & - 'sp median(d2, 2), odd') - -end subroutine - - - - -end program diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp new file mode 100644 index 000000000..3eaca6e00 --- /dev/null +++ b/src/tests/stats/test_median.fypp @@ -0,0 +1,196 @@ +#:include "common.fypp" +#:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES +program test_median + use stdlib_error, only: check + use stdlib_kinds, only: sp, dp, qp, & + int8, int16, int32, int64 + use stdlib_stats, only: median, mean + use,intrinsic :: ieee_arithmetic, only : ieee_is_nan + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 2000 * epsilon(1._dp) + real(qp), parameter :: qptol = 2000 * epsilon(1._qp) + + #:for k1,t1 in IR_KINDS_TYPES + call test_median_${k1}$() + #:endfor + +contains + +#:for k1,t1 in INT_KINDS_TYPES +subroutine test_median_${k1}$() + ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) + + allocate(d0(0)) + call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0') + + d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] + call check(mod(size(d1), 2) == 0, 'test_median_${k1}$: should be an even number') + + d2 = reshape(d1, [3, 4]) + call check(mod(size(d2), 2) == 0, 'test_median_${k1}$: should be an even number') + + d3 = reshape(d1, [2, 3, 2]) + call check(mod(size(d3), 2) == 0, 'test_median_${k1}$: should be an even number') + + !median_all + call check( ieee_is_nan(median(d0)), '${k1}$ median(d0)' ) + call check( ieee_is_nan(median(d0, .false.)), '${k1}$ median(d0, .false.)' ) + call check( ieee_is_nan(median(d1, .false.)), '${k1}$ median(d1, .false.)' ) + call check( ieee_is_nan(median(d2, .false.)), '${k1}$ median(d2, .false.)' ) + call check( ieee_is_nan(median(d3, .false.)), '${k1}$ median(d3, .false.)' ) + + call check( abs(median(d1) - 1.5_dp) < dptol, '${k1}$ median(d1), even') + call check( abs(median(d2) - 1.5_dp) < dptol, '${k1}$ median(d2), even') + call check( abs(median(d3) - 1.5_dp) < dptol, '${k1}$ median(d3), even') + + !median + call check( ieee_is_nan(median(d0, 1)), '${k1}$ median(d0, 1)' ) + call check( ieee_is_nan(median(d0, 1, .false.)), '${k1}$ median(d0, 1, .false.)' ) + call check( ieee_is_nan(median(d1, 1, .false.)), '${k1}$ median(d1, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 1, .false.))), '${k1}$ median(d2, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 2, .false.))), '${k1}$ median(d2, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 1, .false.))), '${k1}$ median(d3, 1, .false.)' ) + call check( any(ieee_is_nan(median(d3, 2, .false.))), '${k1}$ median(d3, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 3, .false.))), '${k1}$ median(d3, 3, .false.)' ) + + call check( abs(median(d1, 1) - 1.5_dp) < dptol, '${k1}$ median(d1, 1), even') + call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp])) < dptol, & + '${k1}$ median(d2, 1)') + call check( sum(abs(median(d2, 2) - [3.5_dp, 1.5_dp, 3._dp])) < dptol, & + '${k1}$ median(d2, 2)') + + !median_mask_all + call check( ieee_is_nan(median(d0, d0 > 0)), '${k1}$ median(d0, d0 > 0)' ) + call check( ieee_is_nan(median(d1, d1 > huge(d1))), '${k1}$ median(d1, d1 > huge(d1))' ) + call check( ieee_is_nan(median(d2, d2 > huge(d2))), '${k1}$ median(d2, d2 > huge(d2))' ) + call check( ieee_is_nan(median(d3, d3 > huge(d3))), '${k1}$ median(d3, d3 > huge(d3))' ) + call check( (median(d1, d1 > 0) - 7._dp) < dptol, '${k1}$ median(d1, d1 > 0)' ) + call check( (median(d2, d2 > 0) - 7._dp) < dptol, '${k1}$ median(d2, d2 > 0)' ) + call check( (median(d3, d3 > 0) - 7._dp) < dptol, '${k1}$ median(d3, d3 > 0)' ) + + !median mask + call check( ieee_is_nan(median(d0, 1, d0 > 0)), '${k1}$ median(d0, 1, d0 > 0)' ) + call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), '${k1}$ median(d1, 1, d1 > huge(d1))' ) + call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), '${k1}$ median(d2, 1, d2 > huge(d2))' ) + call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), '${k1}$ median(d3, 1, d3 > huge(d3))' ) + call check( (median(d1, 1, d1 > 0) - 7._dp) < dptol, '${k1}$ median(d1, 1, d1 >0)') + call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._dp, 6._dp, 8._dp, 10.5_dp] ) )) & + < dptol, '${k1}$ median(d2, 1, d2 > 0 )') + call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_dp, 2._dp, 14.5_dp] ) )) & + < dptol, '${k1}$ median(d2, 2, d2 > 0 )') + call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), '${k1}$ median(d3, 1, d3 > 0)') + + !odd number + d1 = [d1, 20_${k1}$] + call check(mod(size(d1), 2) == 1, 'test_median_${k1}$: should be an odd number') + + call check( abs(median(d1) - 2._dp) < dptol, '${k1}$ median(d1), odd') + + d2 = reshape(d1, [3, 5], pad = [0_${k1}$]) + call check(mod(size(d2), 2) == 1, 'test_median_${k1}$: should be an odd number') + call check( abs(median(d2) - 1._dp) < dptol, '${k1}$ median(d2), odd') + + d3 = reshape(d1, [1, 3, 5], pad = [0_${k1}$]) + call check(mod(size(d3), 2) == 1, 'test_median_${k1}$: should be an odd number') + call check( abs(median(d3) - 1._dp) < dptol, '${k1}$ median(d3), odd') + + call check( abs(median(d1, 1) - 2._dp) < dptol, '${k1}$ median(d1, 1)') + call check( sum(abs(median(d2, 1) - [2._dp, -4._dp, 7._dp, 1._dp, 0._dp])) < dptol, & + '${k1}$ median(d2, 1), odd') + call check( sum(abs(median(d2, 2) - [7._dp, 1._dp, 0._dp])) < dptol, & + '${k1}$ median(d2, 2), odd') + +end subroutine +#:endfor + +#:for k1,t1 in REAL_KINDS_TYPES +subroutine test_median_${k1}$() + ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) + + allocate(d0(0)) + call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0') + + d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] + call check(mod(size(d1), 2) == 0, 'stest_median_${k1}$: hould be an even number') + + d2 = reshape(d1, [3, 4]) + call check(mod(size(d2), 2) == 0, 'test_median_${k1}$: should be an even number') + + d3 = reshape(d1, [2, 3, 2]) + call check(mod(size(d3), 2) == 0, 'test_median_${k1}$: should be an even number') + + !median_all + call check( ieee_is_nan(median(d0)), '${k1}$ median(d0)' ) + call check( ieee_is_nan(median(d0, .false.)), '${k1}$ median(d0, .false.)' ) + call check( ieee_is_nan(median(d1, .false.)), '${k1}$ median(d1, .false.)' ) + call check( ieee_is_nan(median(d2, .false.)), '${k1}$ median(d2, .false.)' ) + call check( ieee_is_nan(median(d3, .false.)), '${k1}$ median(d3, .false.)' ) + + call check( abs(median(d1) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d1), even') + call check( abs(median(d2) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d2), even') + call check( abs(median(d3) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d3), even') + + !median + call check( ieee_is_nan(median(d0, 1)), '${k1}$ median(d0, 1)' ) + call check( ieee_is_nan(median(d0, 1, .false.)), '${k1}$ median(d0, 1, .false.)' ) + call check( ieee_is_nan(median(d1, 1, .false.)), '${k1}$ median(d1, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 1, .false.))), '${k1}$ median(d2, 1, .false.)' ) + call check( any(ieee_is_nan(median(d2, 2, .false.))), '${k1}$ median(d2, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 1, .false.))), '${k1}$ median(d3, 1, .false.)' ) + call check( any(ieee_is_nan(median(d3, 2, .false.))), '${k1}$ median(d3, 2, .false.)' ) + call check( any(ieee_is_nan(median(d3, 3, .false.))), '${k1}$ median(d3, 3, .false.)' ) + + call check( abs(median(d1, 1) - 1.5_${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1), even') + call check( sum(abs(median(d2, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$])) < ${k1}$tol, & + '${k1}$ median(d2, 1)') + call check( sum(abs(median(d2, 2) - [3.5_${k1}$, 1.5_${k1}$, 3._${k1}$])) < ${k1}$tol, & + '${k1}$ median(d2, 2)') + + !median_mask_all + call check( ieee_is_nan(median(d0, d0 > 0)), '${k1}$ median(d0, d0 > 0)' ) + call check( ieee_is_nan(median(d1, d1 > huge(d1))), '${k1}$ median(d1, d1 > huge(d1))' ) + call check( ieee_is_nan(median(d2, d2 > huge(d2))), '${k1}$ median(d2, d2 > huge(d2))' ) + call check( ieee_is_nan(median(d3, d3 > huge(d3))), '${k1}$ median(d3, d3 > huge(d3))' ) + call check( (median(d1, d1 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d1, d1 > 0)' ) + call check( (median(d2, d2 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d2, d2 > 0)' ) + call check( (median(d3, d3 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d3, d3 > 0)' ) + + !median mask + call check( ieee_is_nan(median(d0, 1, d0 > 0)), '${k1}$ median(d0, 1, d0 > 0)' ) + call check( ieee_is_nan(median(d1, 1, d1 > huge(d1))), '${k1}$ median(d1, 1, d1 > huge(d1))' ) + call check( any(ieee_is_nan(median(d2, 1, d2 > huge(d2)))), '${k1}$ median(d2, 1, d2 > huge(d2))' ) + call check( any(ieee_is_nan(median(d3, 1, d3 > huge(d3)))), '${k1}$ median(d3, 1, d3 > huge(d3))' ) + call check( (median(d1, 1, d1 > 0) - 7._${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1, d1 >0)') + call check( sum(abs( (median(d2, 1, d2 > 0) - [ 6._${k1}$, 6._${k1}$, 8._${k1}$, 10.5_${k1}$] ) )) & + < ${k1}$tol, '${k1}$ median(d2, 1, d2 > 0 )') + call check( sum(abs( (median(d2, 2, d2 > 0) - [ 8.5_${k1}$, 2._${k1}$, 14.5_${k1}$] ) )) & + < ${k1}$tol, '${k1}$ median(d2, 2, d2 > 0 )') + call check( any(ieee_is_nan(median(d3, 1, d3 > 0))), '${k1}$ median(d3, 1, d3 > 0)') + + !odd number + d1 = [d1, 20._${k1}$] + call check(mod(size(d1), 2) == 1, 'test_median_${k1}$: should be an odd number') + + call check( abs(median(d1) - 2._${k1}$) < ${k1}$tol, '${k1}$ median(d1), odd') + + d2 = reshape(d1, [3, 5], pad = [0._${k1}$]) + call check(mod(size(d2), 2) == 1, 'test_median_${k1}$: should be an odd number') + call check( abs(median(d2) - 1._${k1}$) < ${k1}$tol, '${k1}$ median(d2), odd') + + d3 = reshape(d1, [1, 3, 5], pad = [0._${k1}$]) + call check(mod(size(d3), 2) == 1, 'test_median_${k1}$: should be an odd number') + call check( abs(median(d3) - 1._${k1}$) < ${k1}$tol, '${k1}$ median(d3), odd') + + call check( abs(median(d1, 1) - 2._${k1}$) < ${k1}$tol, '${k1}$ median(d1, 1)') + call check( sum(abs(median(d2, 1) - [2._${k1}$, -4._${k1}$, 7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol, & + '${k1}$ median(d2, 1), odd') + call check( sum(abs(median(d2, 2) - [7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol, & + '${k1}$ median(d2, 2), odd') + +end subroutine +#:endfor + + + +end program From 3342c6a40119c24974b80315458eaa79fb02c298 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 11 Jun 2021 21:31:43 +0200 Subject: [PATCH 12/30] median: add comments to common.fypp --- src/common.fypp | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index cd4864927..80424db9b 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -147,7 +147,20 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endcall #:enddef - +#! +#! Generates an array rank suffix for subarrays along a dimension +#! +#! Args: +#! rank (int): Rank of the original variable +#! dim (int): Dimension of the variable +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! select_subvector(5, 2) +#! -> (j_, :, j___, j____, j_____) +#! #:def select_subvector(rank, idim) #:assert rank > 0 #:call join_lines(joinstr=", ", prefix="(", suffix=")") @@ -161,6 +174,20 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endcall #:enddef +#! +#! Generates an array rank suffix for arrays +#! +#! Args: +#! rank (int): Rank of the original variable +#! idim (int): Dimension of the variable +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! select_subvector(5, 2) +#! -> (j_, j___, j____, j_____) +#! #:def reduce_subvector(rank, idim) #:assert rank > 0 #:if rank > 1 @@ -175,6 +202,4 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endif #:enddef - - #:endmute From 7e3111e13769daff3c6dddb0846536436e563370 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 16 Jun 2021 18:39:02 +0200 Subject: [PATCH 13/30] median: replace sort to ord_sort --- src/stdlib_stats.fypp | 2 +- src/stdlib_stats_median.fypp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 5a3ed278b..7c05ef0db 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -313,7 +313,7 @@ module stdlib_stats #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS #:set RName = rname("median_all",rank, t1, k1, o1) - pure module function ${RName}$ (x, mask) result(res) + module function ${RName}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index d7ec66af9..3ffc57b0e 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -8,7 +8,7 @@ submodule (stdlib_stats) stdlib_stats_median use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan use stdlib_error, only: error_stop use stdlib_optval, only: optval - use stdlib_sorting, only: sort => sort + use stdlib_sorting, only:sort => ord_sort implicit none contains @@ -16,7 +16,7 @@ contains #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS #:set RName = rname("median_all",rank, t1, k1, o1) - pure module function ${RName}$ (x, mask) result(res) + module function ${RName}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res From dfea79d517cc147ee8ff545b661a8c1effe61eb4 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 17:10:42 +0200 Subject: [PATCH 14/30] update specs --- doc/specs/stdlib_stats.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index 944a111d4..3517180e4 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -189,14 +189,14 @@ Experimental Returns the median of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`. After that the elements are sorted in an increasing order, e.g. `array_sorted = -sort(array)`, the median of the elements of `array` are defined as, if -`n = size(array)` is an even number: +sort(array)`, the median of the elements of `array` is defined as the "middle" +element. If `n = size(array)` is an even number, the median is: ``` median(array) = array_sorted( floor( (n + 1) / 2.)) ``` -or as if `n` is an odd number: +and if `n` is an odd number, the median is: ``` median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) ) @@ -217,7 +217,7 @@ Generic subroutine ### Arguments -`array`: Shall be an array of type `integer` or `real`. +`array`: Shall be an array of type `integer` or `real`. `dim`: Shall be a scalar of type `integer` with a value in the range from 1 to `n`, where `n` is the rank of `array`. @@ -225,7 +225,7 @@ Generic subroutine ### Return value -If `array` is of type `real`, the result is of the same type as `array`. +If `array` is of type `real`, the result is of type `real` with the same kind as `array`. If `array` is of type `integer`, the result is of type `real(dp)`. If `dim` is absent, a scalar with the median of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned. From 8f4a57f0a9ab3f3b9b1f0a775fd1ca22b9775290 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 17:30:58 +0200 Subject: [PATCH 15/30] modif reduce_subvector + select_subvector --- src/common.fypp | 18 ++++++++++-------- src/stdlib_stats_median.fypp | 16 ++++++++-------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index 9f56cad8e..2516aaf79 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -160,6 +160,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! Generates an array rank suffix for subarrays along a dimension #! #! Args: +#! varname (str): Name of the variable to be used as origin #! rank (int): Rank of the original variable #! dim (int): Dimension of the variable #! @@ -167,18 +168,18 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! Array rank suffix string enclosed in braces #! #! E.g., -#! select_subvector(5, 2) +#! select_subvector('j', 5, 2) #! -> (j_, :, j___, j____, j_____) #! -#:def select_subvector(rank, idim) +#:def select_subvector(varname, rank, idim) #:assert rank > 0 #:call join_lines(joinstr=", ", prefix="(", suffix=")") #:for i in range(1, idim) - ${"j" + "_" * (i) }$ + ${varname}$${ "_" * (i) }$ #:endfor : #:for i in range(idim + 1, rank + 1) - ${"j" + "_" * (i) }$ + ${varname}$${ "_" * (i) }$ #:endfor #:endcall #:enddef @@ -187,6 +188,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! Generates an array rank suffix for arrays #! #! Args: +#! varname (str): Name of the variable to be used as origin #! rank (int): Rank of the original variable #! idim (int): Dimension of the variable #! @@ -194,18 +196,18 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! Array rank suffix string enclosed in braces #! #! E.g., -#! select_subvector(5, 2) +#! reduce_subvector('j', 5, 2) #! -> (j_, j___, j____, j_____) #! -#:def reduce_subvector(rank, idim) +#:def reduce_subvector(varname, rank, idim) #:assert rank > 0 #:if rank > 1 #:call join_lines(joinstr=", ", prefix="(", suffix=")") #:for i in range(1, idim) - ${"j" + "_" * (i) }$ + ${varname}$${ "_" * (i) }$ #:endfor #:for i in range(idim + 1, rank + 1) - ${"j" + "_" * (i) }$ + ${varname}$${ "_" * (i) }$ #:endfor #:endcall #:endif diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 3ffc57b0e..b4a9ade7d 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -86,18 +86,18 @@ contains #:for fj in range(fi+1, rank+1) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - x_tmp(:) = x${select_subvector(rank, fi)}$ + x_tmp(:) = x${select_subvector('j', rank, fi)}$ call sort(x_tmp) if (mod(n, 2) == 0) then - res${reduce_subvector(rank, fi)}$ = & + res${reduce_subvector('j', rank, fi)}$ = & #:if t1[0] == 'r' sum(x_tmp(c:c+1)) / 2._${o1}$ #:else sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$ #:endif else - res${reduce_subvector(rank, fi)}$ = x_tmp(c) + res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) end if #:for fj in range(1, rank) end do @@ -180,25 +180,25 @@ contains #:for fj in range(fi+1, rank+1) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - x_tmp = pack(x${select_subvector(rank, fi)}$, & - mask${select_subvector(rank, fi)}$) + x_tmp = pack(x${select_subvector('j', rank, fi)}$, & + mask${select_subvector('j', rank, fi)}$) call sort(x_tmp) n = size(x_tmp, kind=int64) c = floor( (n + 1) / 2._${o1}$, kind=int64 ) if (n == 0) then - res${reduce_subvector(rank, fi)}$ = & + res${reduce_subvector('j', rank, fi)}$ = & ieee_value(1._${o1}$, ieee_quiet_nan) else if (mod(n, 2_int64) == 0) then - res${reduce_subvector(rank, fi)}$ = & + res${reduce_subvector('j', rank, fi)}$ = & #:if t1[0] == 'r' sum(x_tmp(c:c+1)) / 2._${o1}$ #:else sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$ #:endif else if (mod(n, 2_int64) == 1) then - res${reduce_subvector(rank, fi)}$ = x_tmp(c) + res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) end if deallocate(x_tmp) From 4f43ce5a97354340c7132ac84f22faaa9a7aa3b4 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 17:36:13 +0200 Subject: [PATCH 16/30] add comment regarding ord_sort --- doc/specs/stdlib_stats.md | 2 +- src/stdlib_stats_median.fypp | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index 3517180e4..016a2a1d5 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -202,7 +202,7 @@ and if `n` is an odd number, the median is: median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) ) ``` -The array is sorted using the subroutine `sort` provided by the `stdlib_sorting` +The array is sorted using the subroutine `ord_sort` provided by the `stdlib_sorting` module. ### Syntax diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index b4a9ade7d..5c899256d 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -8,6 +8,9 @@ submodule (stdlib_stats) stdlib_stats_median use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan use stdlib_error, only: error_stop use stdlib_optval, only: optval + ! Use "ord_sort" rather than "sort" because the former can be much faster for arrays + ! that are already partly sorted. While it is slightly slower for random arrays, + ! ord_sort seems a better overall choice. use stdlib_sorting, only:sort => ord_sort implicit none From 6644140788e7f7cfb8c7caeb3c7b1787052f7f02 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 18:16:25 +0200 Subject: [PATCH 17/30] add comment + combined fypp loops --- src/stdlib_stats_median.fypp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 5c899256d..85d91b7dd 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -83,12 +83,10 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - #:for fj in range(1, fi) + ! Loop over every dimension of the array except "dim" + #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - #:for fj in range(fi+1, rank+1) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor x_tmp(:) = x${select_subvector('j', rank, fi)}$ call sort(x_tmp) @@ -177,12 +175,10 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - #:for fj in range(1, fi) + ! Loop over every dimension of the array except "dim" + #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - #:for fj in range(fi+1, rank+1) - do j${"_" * fj}$ = 1, size(x, ${fj}$) - #:endfor x_tmp = pack(x${select_subvector('j', rank, fi)}$, & mask${select_subvector('j', rank, fi)}$) call sort(x_tmp) From 2ca833a392b16c364f569f94d0747baa19670083 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 18:35:50 +0200 Subject: [PATCH 18/30] update hyperlinks --- src/stdlib_stats.fypp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 7c05ef0db..c7d376e28 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -19,7 +19,7 @@ module stdlib_stats !! version: experimental !! !! Pearson correlation of array elements - !! ([Specification](../page/specs/stdlib_stats.html#description)) + !! ([Specification](../page/specs/stdlib_stats.html#corr-pearson-correlation-of-array-elements)) #:for k1, t1 in RC_KINDS_TYPES #:set RName = rname("corr",1, t1, k1) module function ${RName}$(x, dim, mask) result(res) @@ -111,7 +111,7 @@ module stdlib_stats !! version: experimental !! !! Covariance of array elements - !! ([Specification](../page/specs/stdlib_stats.html#description_1)) + !! ([Specification](../page/specs/stdlib_stats.html#cov-covariance-of-array-elements)) #:for k1, t1 in RC_KINDS_TYPES #:set RName = rname("cov",1, t1, k1) module function ${RName}$(x, dim, mask, corrected) result(res) @@ -210,7 +210,7 @@ module stdlib_stats !! version: experimental !! !! Mean of array elements - !! ([Specification](../page/specs/stdlib_stats.html#description_2)) + !! ([Specification](../page/specs/stdlib_stats.html#mean-mean-of-array-elements)) #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("mean_all",rank, t1, k1) @@ -305,11 +305,12 @@ module stdlib_stats end interface mean + interface median !! version: experimental !! !! Median of array elements - !! ([Specification](../page/specs/stdlib_stats.html#median)) + !! ([Specification](../page/specs/stdlib_stats.html#median-median-of-array-elements)) #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS #:set RName = rname("median_all",rank, t1, k1, o1) @@ -363,7 +364,7 @@ module stdlib_stats !! version: experimental !! !! Variance of array elements - !! ([Specification](../page/specs/stdlib_stats.html#description_4)) + !! ([Specification](../page/specs/stdlib_stats.html#var-variance-of-array-elements)) #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS @@ -472,7 +473,7 @@ module stdlib_stats !! version: experimental !! !! Central moment of array elements - !! ([Specification](../page/specs/stdlib_stats.html#description_3)) + !! ([Specification](../page/specs/stdlib_stats.html#moment-central-moments-of-array-elements)) #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_all",rank, t1, k1) From e0f68ed61f44893fe2874ca11fb3ba713fa6338f Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 18:40:36 +0200 Subject: [PATCH 19/30] add link in stdlib_stats.md --- doc/specs/stdlib_stats.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index 016a2a1d5..45b2db1ce 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -202,8 +202,8 @@ and if `n` is an odd number, the median is: median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) ) ``` -The array is sorted using the subroutine `ord_sort` provided by the `stdlib_sorting` -module. +The array is sorted using the subroutine `[[stdlib_sorting(module):ord_sort(interface)]]` +provided by the `[[stdlib_sorting(module)]]` module. ### Syntax From 3eb4e4afadb6c2751475b92a82cd8180e1d10ac9 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sun, 4 Jul 2021 19:15:36 +0200 Subject: [PATCH 20/30] formatting --- src/stdlib_stats_median.fypp | 70 ++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 85d91b7dd..99ffaaff0 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -87,19 +87,19 @@ contains #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - x_tmp(:) = x${select_subvector('j', rank, fi)}$ - call sort(x_tmp) - - if (mod(n, 2) == 0) then - res${reduce_subvector('j', rank, fi)}$ = & - #:if t1[0] == 'r' - sum(x_tmp(c:c+1)) / 2._${o1}$ - #:else - sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$ - #:endif - else - res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) - end if + x_tmp(:) = x${select_subvector('j', rank, fi)}$ + call sort(x_tmp) + + if (mod(n, 2) == 0) then + res${reduce_subvector('j', rank, fi)}$ = & + #:if t1[0] == 'r' + sum(x_tmp(c:c+1)) / 2._${o1}$ + #:else + sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$ + #:endif + else + res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) + end if #:for fj in range(1, rank) end do #:endfor @@ -179,28 +179,28 @@ contains #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) do j${"_" * fj}$ = 1, size(x, ${fj}$) #:endfor - x_tmp = pack(x${select_subvector('j', rank, fi)}$, & - mask${select_subvector('j', rank, fi)}$) - call sort(x_tmp) - - n = size(x_tmp, kind=int64) - c = floor( (n + 1) / 2._${o1}$, kind=int64 ) - - if (n == 0) then - res${reduce_subvector('j', rank, fi)}$ = & - ieee_value(1._${o1}$, ieee_quiet_nan) - else if (mod(n, 2_int64) == 0) then - res${reduce_subvector('j', rank, fi)}$ = & - #:if t1[0] == 'r' - sum(x_tmp(c:c+1)) / 2._${o1}$ - #:else - sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$ - #:endif - else if (mod(n, 2_int64) == 1) then - res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) - end if - - deallocate(x_tmp) + x_tmp = pack(x${select_subvector('j', rank, fi)}$, & + mask${select_subvector('j', rank, fi)}$) + call sort(x_tmp) + + n = size(x_tmp, kind=int64) + c = floor( (n + 1) / 2._${o1}$, kind=int64 ) + + if (n == 0) then + res${reduce_subvector('j', rank, fi)}$ = & + ieee_value(1._${o1}$, ieee_quiet_nan) + else if (mod(n, 2_int64) == 0) then + res${reduce_subvector('j', rank, fi)}$ = & + #:if t1[0] == 'r' + sum(x_tmp(c:c+1)) / 2._${o1}$ + #:else + sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$ + #:endif + else if (mod(n, 2_int64) == 1) then + res${reduce_subvector('j', rank, fi)}$ = x_tmp(c) + end if + + deallocate(x_tmp) #:for fj in range(1, rank) end do #:endfor From de2c41918f1f5d3e01680a9892b99df1ad6a4a39 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 04:04:01 -0400 Subject: [PATCH 21/30] Apply suggestions from code review Co-authored-by: Ivan Pribec --- doc/specs/stdlib_stats.md | 6 +++--- src/stdlib_stats_median.fypp | 2 +- src/tests/stats/test_median.fypp | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index 45b2db1ce..e6947f425 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -188,9 +188,9 @@ Experimental Returns the median of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`. -After that the elements are sorted in an increasing order, e.g. `array_sorted = -sort(array)`, the median of the elements of `array` is defined as the "middle" -element. If `n = size(array)` is an even number, the median is: +The median of the elements of `array` is defined as the "middle" +element, after that the elements are sorted in an increasing order, e.g. `array_sorted = +sort(array)`. If `n = size(array)` is an even number, the median is: ``` median(array) = array_sorted( floor( (n + 1) / 2.)) diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 99ffaaff0..3d78385ad 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -11,7 +11,7 @@ submodule (stdlib_stats) stdlib_stats_median ! Use "ord_sort" rather than "sort" because the former can be much faster for arrays ! that are already partly sorted. While it is slightly slower for random arrays, ! ord_sort seems a better overall choice. - use stdlib_sorting, only:sort => ord_sort + use stdlib_sorting, only: sort => ord_sort implicit none contains diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp index 3eaca6e00..9f4d70c4b 100644 --- a/src/tests/stats/test_median.fypp +++ b/src/tests/stats/test_median.fypp @@ -5,7 +5,7 @@ program test_median use stdlib_kinds, only: sp, dp, qp, & int8, int16, int32, int64 use stdlib_stats, only: median, mean - use,intrinsic :: ieee_arithmetic, only : ieee_is_nan + use, intrinsic :: ieee_arithmetic, only : ieee_is_nan implicit none real(sp), parameter :: sptol = 1000 * epsilon(1._sp) real(dp), parameter :: dptol = 2000 * epsilon(1._dp) @@ -112,7 +112,7 @@ subroutine test_median_${k1}$() call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0') d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] - call check(mod(size(d1), 2) == 0, 'stest_median_${k1}$: hould be an even number') + call check(mod(size(d1), 2) == 0, 'test_median_${k1}$: should be an even number') d2 = reshape(d1, [3, 4]) call check(mod(size(d2), 2) == 0, 'test_median_${k1}$: should be an even number') From 9bbcb7482e88ea13a2b4b43f58a3bcf68b83acb6 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 18:18:54 +0200 Subject: [PATCH 22/30] median: reorder fypp variable --- src/stdlib_stats.fypp | 2 +- src/stdlib_stats_median.fypp | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index c7d376e28..8b9b2dd63 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -2,7 +2,7 @@ #:set RANKS = range(1, MAXRANK + 1) #:set REDRANKS = range(2, MAXRANK + 1) #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -#:set IR_KINDS_TYPES_OUTPUT = list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) + list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) +#:set IR_KINDS_TYPES_OUTPUT = list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) + list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) module stdlib_stats !! Provides support for various statistical methods. This includes currently !! descriptive statistics diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index 3d78385ad..efe5ebdd5 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -1,7 +1,6 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) -#:set IR_KINDS_TYPES_OUTPUT = list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) + list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) - +#:set IR_KINDS_TYPES_OUTPUT = list(zip(INT_KINDS,INT_TYPES, ['dp']*len(INT_KINDS))) + list(zip(REAL_KINDS, REAL_TYPES, REAL_KINDS)) submodule (stdlib_stats) stdlib_stats_median From f17b890e7d42ff192ce7af8d1730de6b672b7ed7 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 18:50:11 +0200 Subject: [PATCH 23/30] median: replace _ by numbers --- src/common.fypp | 12 ++++++------ src/stdlib_stats_median.fypp | 8 ++++---- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index 2516aaf79..bfe7743ec 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -169,17 +169,17 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! #! E.g., #! select_subvector('j', 5, 2) -#! -> (j_, :, j___, j____, j_____) +#! -> (j1, :, j3, j4, j5) #! #:def select_subvector(varname, rank, idim) #:assert rank > 0 #:call join_lines(joinstr=", ", prefix="(", suffix=")") #:for i in range(1, idim) - ${varname}$${ "_" * (i) }$ + ${varname}$${i}$ #:endfor : #:for i in range(idim + 1, rank + 1) - ${varname}$${ "_" * (i) }$ + ${varname}$${i}$ #:endfor #:endcall #:enddef @@ -197,17 +197,17 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! #! E.g., #! reduce_subvector('j', 5, 2) -#! -> (j_, j___, j____, j_____) +#! -> (j1, j3, j4, j5) #! #:def reduce_subvector(varname, rank, idim) #:assert rank > 0 #:if rank > 1 #:call join_lines(joinstr=", ", prefix="(", suffix=")") #:for i in range(1, idim) - ${varname}$${ "_" * (i) }$ + ${varname}$${i}$ #:endfor #:for i in range(idim + 1, rank + 1) - ${varname}$${ "_" * (i) }$ + ${varname}$${i}$ #:endfor #:endcall #:endif diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index efe5ebdd5..dcdb1c814 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -64,7 +64,7 @@ contains integer :: c, n #:if rank > 1 #:for fj in range(1, rank+1) - integer :: j${"_" * fj}$ + integer :: j${fj}$ #:endfor #:endif ${t1}$, allocatable :: x_tmp(:) @@ -84,7 +84,7 @@ contains case(${fi}$) ! Loop over every dimension of the array except "dim" #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) - do j${"_" * fj}$ = 1, size(x, ${fj}$) + do j${fj}$ = 1, size(x, ${fj}$) #:endfor x_tmp(:) = x${select_subvector('j', rank, fi)}$ call sort(x_tmp) @@ -162,7 +162,7 @@ contains integer(kind = int64) :: c, n #:if rank > 1 #:for fj in range(1, rank+1) - integer :: j${"_" * fj}$ + integer :: j${fj}$ #:endfor #:endif ${t1}$, allocatable :: x_tmp(:) @@ -176,7 +176,7 @@ contains case(${fi}$) ! Loop over every dimension of the array except "dim" #:for fj in list(range(1, fi)) + list(range(fi+1, rank+1)) - do j${"_" * fj}$ = 1, size(x, ${fj}$) + do j${fj}$ = 1, size(x, ${fj}$) #:endfor x_tmp = pack(x${select_subvector('j', rank, fi)}$, & mask${select_subvector('j', rank, fi)}$) From 391c658858051af353d7425c53b5c2b60d8dc296 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 18:51:30 +0200 Subject: [PATCH 24/30] median: add in common.fypp where it is used for median case --- src/common.fypp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/common.fypp b/src/common.fypp index bfe7743ec..360b2da3c 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -171,6 +171,9 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! select_subvector('j', 5, 2) #! -> (j1, :, j3, j4, j5) #! +#! Used, e.g., in +#! stdlib_stats_median.fypp +#! #:def select_subvector(varname, rank, idim) #:assert rank > 0 #:call join_lines(joinstr=", ", prefix="(", suffix=")") @@ -199,6 +202,9 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! reduce_subvector('j', 5, 2) #! -> (j1, j3, j4, j5) #! +#! Used, e.g., in +#! stdlib_stats_median.fypp +#! #:def reduce_subvector(varname, rank, idim) #:assert rank > 0 #:if rank > 1 From afc92a22ac1139f74980428c507b8a64d3afe14d Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 19:20:47 +0200 Subject: [PATCH 25/30] median: add comment in test median --- src/tests/stats/test_median.fypp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp index 9f4d70c4b..f82732d9e 100644 --- a/src/tests/stats/test_median.fypp +++ b/src/tests/stats/test_median.fypp @@ -22,6 +22,7 @@ subroutine test_median_${k1}$() ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) allocate(d0(0)) + !check just to be sure that the setup of d0 is correct call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0') d1 = [10, 2, -3, -4, 6, -6, 7, -8, 9, 0, 1, 20] @@ -109,6 +110,7 @@ subroutine test_median_${k1}$() ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) allocate(d0(0)) + !check just to be sure that the setup of d0 is correct call check(size(d0) == 0, 'test_median_${k1}$: should be of size 0') d1 = [10., 2., -3., -4., 6., -6., 7., -8., 9., 0., 1., 20.] From 4d328dc56bda6d8b58f2360b7fb680ddfb38b4c3 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 19:34:32 +0200 Subject: [PATCH 26/30] median: add warning about naive implementation --- doc/specs/stdlib_stats.md | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index e6947f425..1c0d42944 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -186,7 +186,9 @@ Experimental ### Description -Returns the median of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`. +Returns the median of all the elements of `array`, or of the elements of `array` +along dimension `dim` if provided, and if the corresponding element in `mask` is +`true`. The median of the elements of `array` is defined as the "middle" element, after that the elements are sorted in an increasing order, e.g. `array_sorted = @@ -202,7 +204,8 @@ and if `n` is an odd number, the median is: median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) ) ``` -The array is sorted using the subroutine `[[stdlib_sorting(module):ord_sort(interface)]]` +The current implementation is a quite naive implementation that relies on sorting +the whole array, using the subroutine `[[stdlib_sorting(module):ord_sort(interface)]]` provided by the `[[stdlib_sorting(module)]]` module. ### Syntax From bdb47b790a839a830df901b2f0cff7d8565558ba Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Wed, 21 Jul 2021 23:07:11 +0200 Subject: [PATCH 27/30] median: return NaN when real array contain NaN --- doc/specs/stdlib_stats.md | 2 ++ src/stdlib_stats_median.fypp | 46 +++++++++++++++++++++++++++++--- src/tests/stats/test_median.fypp | 38 +++++++++++++++++++++++++- 3 files changed, 82 insertions(+), 4 deletions(-) diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md index 1c0d42944..b8a4f7f78 100644 --- a/doc/specs/stdlib_stats.md +++ b/doc/specs/stdlib_stats.md @@ -229,12 +229,14 @@ Generic subroutine ### Return value If `array` is of type `real`, the result is of type `real` with the same kind as `array`. +If `array` is of type `real` and contains IEEE `NaN`, the result is IEEE `NaN`. If `array` is of type `integer`, the result is of type `real(dp)`. If `dim` is absent, a scalar with the median of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned. If `mask` is specified, the result is the median of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`. + ### Example ```fortran diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index dcdb1c814..b91ff54cf 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -4,7 +4,7 @@ submodule (stdlib_stats) stdlib_stats_median - use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan use stdlib_error, only: error_stop use stdlib_optval, only: optval ! Use "ord_sort" rather than "sort" because the former can be much faster for arrays @@ -31,6 +31,13 @@ contains return end if + #:if t1[0] == 'r' + if (any(ieee_is_nan(x))) then + res = ieee_value(1._${o1}$, ieee_quiet_nan) + return + end if + #:endif + n = size(x, kind=int64) c = floor( (n + 1) / 2._${o1}$, kind=int64 ) @@ -87,6 +94,19 @@ contains do j${fj}$ = 1, size(x, ${fj}$) #:endfor x_tmp(:) = x${select_subvector('j', rank, fi)}$ + + #:if t1[0] == 'r' + if (any(ieee_is_nan(x_tmp))) then + res${reduce_subvector('j', rank, fi)}$ = & + ieee_value(1._${o1}$, ieee_quiet_nan) + #:if fi == 1 + return + #:else + cycle + #:endif + end if + #:endif + call sort(x_tmp) if (mod(n, 2) == 0) then @@ -127,12 +147,19 @@ contains call error_stop("ERROR (median): shapes of x and mask are different") end if + #:if t1[0] == 'r' + if (any(ieee_is_nan(x))) then + res = ieee_value(1._${o1}$, ieee_quiet_nan) + return + end if + #:endif + x_tmp = pack(x, mask) call sort(x_tmp) - n = size(x_tmp, kind=int64 ) - c = floor( (n + 1) / 2._${o1}$, kind=int64 ) + n = size(x_tmp, kind=int64) + c = floor( (n + 1) / 2._${o1}$, kind=int64) if (n == 0) then res = ieee_value(1._${o1}$, ieee_quiet_nan) @@ -180,6 +207,19 @@ contains #:endfor x_tmp = pack(x${select_subvector('j', rank, fi)}$, & mask${select_subvector('j', rank, fi)}$) + + #:if t1[0] == 'r' + if (any(ieee_is_nan(x_tmp))) then + res${reduce_subvector('j', rank, fi)}$ = & + ieee_value(1._${o1}$, ieee_quiet_nan) + #:if rank == 1 + return + #:else + cycle + #:endif + end if + #:endif + call sort(x_tmp) n = size(x_tmp, kind=int64) diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp index f82732d9e..09962da27 100644 --- a/src/tests/stats/test_median.fypp +++ b/src/tests/stats/test_median.fypp @@ -5,7 +5,7 @@ program test_median use stdlib_kinds, only: sp, dp, qp, & int8, int16, int32, int64 use stdlib_stats, only: median, mean - use, intrinsic :: ieee_arithmetic, only : ieee_is_nan + use, intrinsic :: ieee_arithmetic, only : ieee_is_nan, ieee_value, ieee_quiet_nan implicit none real(sp), parameter :: sptol = 1000 * epsilon(1._sp) real(dp), parameter :: dptol = 2000 * epsilon(1._dp) @@ -108,6 +108,7 @@ end subroutine #:for k1,t1 in REAL_KINDS_TYPES subroutine test_median_${k1}$() ${t1}$, allocatable :: d0(:), d1(:), d2(:,:), d3(:,:,:) + ${t1}$, allocatable :: tmp1(:) allocate(d0(0)) !check just to be sure that the setup of d0 is correct @@ -190,6 +191,41 @@ subroutine test_median_${k1}$() call check( sum(abs(median(d2, 2) - [7._${k1}$, 1._${k1}$, 0._${k1}$])) < ${k1}$tol, & '${k1}$ median(d2, 2), odd') + !check IEEE NaN values in array + d1(1) = ieee_value(1._${k1}$, ieee_quiet_nan) + d2(1, 1) = ieee_value(1._${k1}$, ieee_quiet_nan) + d3(1, 1, 1) = ieee_value(1._${k1}$, ieee_quiet_nan) + + !median_all + call check( ieee_is_nan(median(d1)), '${k1}$ median(d1), should be NaN') + call check( ieee_is_nan(median(d2)), '${k1}$ median(d2), should be NaN') + call check( ieee_is_nan(median(d3)), '${k1}$ median(d3), should be NaN') + + !median + call check( any(ieee_is_nan(median(d2, 1))), '${k1}$ median(d2, 1) should contain at least 1 NaN' ) + call check( any(ieee_is_nan(median(d2, 2))), '${k1}$ median(d2, 2) should contain at least 1 NaN' ) + + call check( any(ieee_is_nan(median(d3, 1))), '${k1}$ median(d3, 1) should contain at least 1 NaN' ) + call check( any(ieee_is_nan(median(d3, 2))), '${k1}$ median(d3, 2) should contain at least 1 NaN' ) + call check( any(ieee_is_nan(median(d3, 3))), '${k1}$ median(d3, 3) should contain at least 1 NaN' ) + + !median_mask_all + call check( ieee_is_nan(median(d1, d1 > 0)), '${k1}$ median(d1, d1 > 0) should be NaN' ) + call check( ieee_is_nan(median(d2, d2 > 0)), '${k1}$ median(d2, d2 > 0) should be NaN' ) + call check( ieee_is_nan(median(d3, d3 > 0)), '${k1}$ median(d3, d3 > 0) should be NaN' ) + + !median mask + call check( ieee_is_nan(median(d1, 1, d1 .ne. 0)), '${k1}$ median(d1, 1, d1.ne.0) should return NaN') + + tmp1 = median(d2, 1, d2 .ne. 0) + call check( any(ieee_is_nan(tmp1)), '${k1}$ median(d2, 1, d2 .ne. 0 ) should contain at least 1 NaN') + call check( tmp1(2) == -4._${k1}$, '${k1}$ tmp1(2) == -4') + call check( any(ieee_is_nan(median(d2, 2, d2 .ne. 0))), & + '${k1}$ median(d2, 2, d2 .ne. 0 ) should contain at least 1 NaN') + + call check( any(ieee_is_nan(median(d3, 1, d3 .ne. 0))), & + '${k1}$ median(d3, 1, d3 .ne. 0 ) should contain at least 1 NaN') + end subroutine #:endfor From 227f021011828d43ea7375b8c4bcaa03b691860c Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 23 Jul 2021 03:11:47 -0400 Subject: [PATCH 28/30] @ivan-pi suggestions from code review Co-authored-by: Ivan Pribec --- src/common.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index 360b2da3c..0c6e86840 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -192,8 +192,8 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #! #! Args: #! varname (str): Name of the variable to be used as origin -#! rank (int): Rank of the original variable -#! idim (int): Dimension of the variable +#! rank (int): Rank of the original array variable +#! idim (int): Dimension of the variable dropped #! #! Returns: #! Array rank suffix string enclosed in braces From 9d38d9d3ed44d51aacabc26258e14afb80c509b9 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 23 Jul 2021 09:44:22 +0200 Subject: [PATCH 29/30] median: rename fypp RName by name --- src/stdlib_stats.fypp | 24 ++++++++++++------------ src/stdlib_stats_median.fypp | 24 ++++++++++++------------ 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index 8b9b2dd63..d52be417b 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -313,47 +313,47 @@ module stdlib_stats !! ([Specification](../page/specs/stdlib_stats.html#median-median-of-array-elements)) #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1, o1) - module function ${RName}$ (x, mask) result(res) + #:set name = rname("median_all",rank, t1, k1, o1) + module function ${name}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1, o1) - module function ${RName}$(x, dim, mask) result(res) + #:set name = rname("median",rank, t1, k1, o1) + module function ${name}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1, o1) - module function ${RName}$(x, mask) result(res) + #:set name = rname('median_mask_all',rank, t1, k1, o1) + module function ${name}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ real(${o1}$) :: res - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1, o1) - module function ${RName}$(x, dim, mask) result(res) + #:set name = rname('median_mask',rank, t1, k1, o1) + module function ${name}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ real(${o1}$) :: res${reduced_shape('x', rank, 'dim')}$ - end function ${RName}$ + end function ${name}$ #:endfor #:endfor diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index b91ff54cf..b35999cf3 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -17,8 +17,8 @@ contains #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median_all",rank, t1, k1, o1) - module function ${RName}$ (x, mask) result(res) + #:set name = rname("median_all",rank, t1, k1, o1) + module function ${name}$ (x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${o1}$) :: res @@ -55,14 +55,14 @@ contains res = x_tmp(c) end if - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname("median",rank, t1, k1, o1) - module function ${RName}$(x, dim, mask) result(res) + #:set name = rname("median",rank, t1, k1, o1) + module function ${name}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask @@ -127,15 +127,15 @@ contains call error_stop("ERROR (median): wrong dimension") end select - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask_all',rank, t1, k1, o1) - module function ${RName}$(x, mask) result(res) + #:set name = rname('median_mask_all',rank, t1, k1, o1) + module function ${name}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ real(${o1}$) :: res @@ -173,14 +173,14 @@ contains res = x_tmp(c) end if - end function ${RName}$ + end function ${name}$ #:endfor #:endfor #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set RName = rname('median_mask',rank, t1, k1, o1) - module function ${RName}$(x, dim, mask) result(res) + #:set name = rname('median_mask',rank, t1, k1, o1) + module function ${name}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ @@ -248,7 +248,7 @@ contains call error_stop("ERROR (median): wrong dimension") end select - end function ${RName}$ + end function ${name}$ #:endfor #:endfor From cbdc4acbad4926f70eea5ed604cf5b64cbfbd801 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 23 Jul 2021 09:50:34 +0200 Subject: [PATCH 30/30] median: replace median_mask_all by median_all_mask --- src/stdlib_stats.fypp | 2 +- src/stdlib_stats_median.fypp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_stats.fypp b/src/stdlib_stats.fypp index d52be417b..038f27752 100644 --- a/src/stdlib_stats.fypp +++ b/src/stdlib_stats.fypp @@ -336,7 +336,7 @@ module stdlib_stats #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set name = rname('median_mask_all',rank, t1, k1, o1) + #:set name = rname('median_all_mask',rank, t1, k1, o1) module function ${name}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp index b35999cf3..5adb7b256 100644 --- a/src/stdlib_stats_median.fypp +++ b/src/stdlib_stats_median.fypp @@ -134,7 +134,7 @@ contains #:for k1, t1, o1 in IR_KINDS_TYPES_OUTPUT #:for rank in RANKS - #:set name = rname('median_mask_all',rank, t1, k1, o1) + #:set name = rname('median_all_mask',rank, t1, k1, o1) module function ${name}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$