From e966e7bacc689356e9cdd971700aadfed334fb1f Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Wed, 5 Feb 2020 23:36:00 +0100 Subject: [PATCH 1/8] addition of variance --- src/CMakeLists.txt | 1 + src/Makefile.manual | 8 +- src/common.fypp | 15 + src/stdlib_experimental_stats.fypp | 96 +++++- src/stdlib_experimental_stats.md | 53 ++++ src/stdlib_experimental_stats_var.fypp | 240 +++++++++++++++ src/tests/stats/CMakeLists.txt | 1 + src/tests/stats/Makefile.manual | 2 +- src/tests/stats/test_var.f90 | 391 +++++++++++++++++++++++++ 9 files changed, 804 insertions(+), 3 deletions(-) create mode 100644 src/stdlib_experimental_stats_var.fypp create mode 100644 src/tests/stats/test_var.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 106fa02ac..df4358220 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ set(fppFiles stdlib_experimental_optval.fypp stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.fypp + stdlib_experimental_stats_var.fypp ) diff --git a/src/Makefile.manual b/src/Makefile.manual index f3772f2bd..7159cfa76 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -5,7 +5,8 @@ SRC = stdlib_experimental_ascii.f90 \ stdlib_experimental_kinds.f90 \ f18estop.f90 \ stdlib_experimental_stats.f90 \ - stdlib_experimental_stats_mean.f90 + stdlib_experimental_stats_mean.f90 \ + stdlib_experimental_stats_var.f90 LIB = libstdlib.a @@ -42,6 +43,11 @@ stdlib_experimental_stats_mean.o: \ stdlib_experimental_optval.o \ stdlib_experimental_kinds.o \ stdlib_experimental_stats.o +stdlib_experimental_stats_var.o: \ + stdlib_experimental_optval.o \ + stdlib_experimental_kinds.o \ + stdlib_experimental_stats.o stdlib_experimental_io.f90: stdlib_experimental_io.fypp stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp +stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp diff --git a/src/common.fypp b/src/common.fypp index 99fed2250..b3e45d3d4 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -113,4 +113,19 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ $:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix) #:enddef +#:def rankindice(varname, varname1, origrank, dim) + #:assert origrank > 0 + #:if origrank > 0 + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, origrank+1) + #:if i == dim + ${varname1}$ + #:else + ${varname}$ + #:endif + #:endfor + #:endcall + #:endif +#:enddef + #:endmute diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index f16fa0a92..ab0582a4d 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -7,7 +7,7 @@ module stdlib_experimental_stats implicit none private ! Public API - public :: mean + public :: mean, var interface mean #:for k1, t1 in RC_KINDS_TYPES @@ -104,4 +104,98 @@ module stdlib_experimental_stats end interface mean + + interface var + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + ${t1}$ :: res + end function var_${rank}$_all_${k1}$_${k1}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_all_${k1}$_dp(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + end function var_${rank}$_all_${k1}$_dp + #:endfor + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_${k1}$_${k1}$(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 var_${rank}$_${k1}$_${k1}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_${k1}$_dp(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 var_${rank}$_${k1}$_dp + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res + end function var_${rank}$_mask_all_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + end function var_${rank}$_mask_all_${k1}$_dp + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_${k1}$_${k1}$(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 var_${rank}$_mask_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_${k1}$_dp(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 var_${rank}$_mask_${k1}$_dp + #:endfor + #:endfor + + end interface var + + end module stdlib_experimental_stats diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index 0c11c6822..00f96d6df 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -3,6 +3,7 @@ ## Implemented * `mean` + * `var` ## `mean` - mean of array elements @@ -47,3 +48,55 @@ program demo_mean reshape(x, [ 2, 3 ] ) > 3.) !returns [ NaN, 4.0, 5.5 ] end program demo_mean ``` + +## `var` - variance of array elements + +### Description + +Returns the variance 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 variance is defined as the best unbiased estimator and is computed as: + +``` + var(x) = 1/(n-1) sum_i (array(i) - mean(array))^2 +``` + +### Syntax + +`result = var(array [, mask])` + +`result = var(array, dim [, mask])` + +### 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 by 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 `double precision`. + +If `dim` is absent, a scalar with the variance 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 `ar +ray` with dimension `dim` dropped is returned. + +If `mask` is specified, the result is the variance 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_var + use stdlib_experimental_stats, only: var + implicit none + real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ] + print *, var(x) !returns 3.5 + print *, var( reshape(x, [ 2, 3 ] )) !returns 3.5 + print *, var( reshape(x, [ 2, 3 ] ), 1) !returns [0.5, 0.5, 0.5] + print *, var( reshape(x, [ 2, 3 ] ), 1,& + reshape(x, [ 2, 3 ] ) > 3.) !returns [0., NaN, 0.5] +end program demo_var +``` + diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp new file mode 100644 index 000000000..675767a0a --- /dev/null +++ b/src/stdlib_experimental_stats_var.fypp @@ -0,0 +1,240 @@ +#:include "common.fypp" + +#:set RANKS = range(1, MAXRANK + 1) + + +submodule (stdlib_experimental_stats) stdlib_experimental_stats_var + + use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use stdlib_experimental_error, only: error_stop + use stdlib_experimental_optval, only: optval + implicit none + +contains + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + ${t1}$ :: res + + ${t1}$ :: n, mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(res, ieee_quiet_nan) + return + end if + + n = real(size(x, kind = int64), ${k1}$) + mean = sum(x) / n + + res = sum((x - mean)**2) / (n - 1._${k1}$) + + end function var_${rank}$_all_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_all_${k1}$_dp(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + + real(dp) :: n, mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(res, ieee_quiet_nan) + return + end if + + n = real(size(x, kind = int64), dp) + mean = sum(real(x, dp)) / n + + res = sum((real(x, dp) - mean)**2) / (n - 1._dp) + + end function var_${rank}$_all_${k1}$_dp + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_${k1}$_${k1}$(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 :: i + ${t1}$ :: n + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(res, ieee_quiet_nan) + return + end if + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(size(x, dim), ${k1}$) + mean = sum(x, dim) / n + do i = 1, size(x, dim) + res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2 + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function var_${rank}$_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_${k1}$_dp(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 :: i + real(dp) :: n + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(res, ieee_quiet_nan) + return + end if + + res = 0._dp + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(size(x, dim), dp) + mean = sum(real(x, dp), dim) / n + do i = 1, size(x, dim) + res = res + (real(x${rankindice(':', 'i', rank, fi )}$, dp) - mean)**2 + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._dp) + + end function var_${rank}$_${k1}$_dp + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res + + ${t1}$ :: n, mean + + n = real(count(mask, kind = int64), ${k1}$) + mean = sum(x, mask) / n + + res = sum((x - mean)**2, mask) / (n - 1._${k1}$) + + end function var_${rank}$_mask_all_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + + real(dp) :: n, mean + + n = real(count(mask, kind = int64), dp) + mean = sum(real(x, dp), mask) / n + + res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp) + + end function var_${rank}$_mask_all_${k1}$_dp + #:endfor + #:endfor + + + #:for k1, t1 in REAL_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_${k1}$_${k1}$(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 :: i + ${t1}$ :: n${reduced_shape('x', rank, 'dim')}$ + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(count(mask, dim), ${k1}$) + mean = sum(x, dim, mask) / n + do i = 1, size(x, dim) + res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& + 0._${k1}$, mask${rankindice(':', 'i', rank, fi)}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function var_${rank}$_mask_${k1}$_${k1}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + module function var_${rank}$_mask_${k1}$_dp(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 :: i + real(dp) :: n${reduced_shape('x', rank, 'dim')}$ + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0._dp + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(count(mask, dim), dp) + mean = sum(real(x, dp), dim, mask) / n + do i = 1, size(x, dim) + res = res + merge((real(x${rankindice(':', 'i', rank, fi )}$, dp) - mean)**2,& + 0._dp, mask${rankindice(':', 'i', rank, fi)}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._dp) + + end function var_${rank}$_mask_${k1}$_dp + #:endfor + #:endfor + +end submodule diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index e86fe23d7..239ac2373 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,4 +1,5 @@ ADDTEST(mean) +ADDTEST(var) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 9faf154cb..1d25a5509 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,3 @@ -PROGS_SRC = test_mean.f90 +PROGS_SRC = test_mean.f90 test_var.f90 include ../Makefile.manual.test.mk diff --git a/src/tests/stats/test_var.f90 b/src/tests/stats/test_var.f90 new file mode 100644 index 000000000..32481615d --- /dev/null +++ b/src/tests/stats/test_var.f90 @@ -0,0 +1,391 @@ +program test_var + use stdlib_experimental_error, only: assert + use stdlib_experimental_kinds, only: sp, dp, int32, int64 + use stdlib_experimental_io, only: loadtxt + use stdlib_experimental_stats, only: mean, var + implicit none + + + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 1000 * epsilon(1._dp) + + integer(int32) :: i321(5) = [1, 2, 3, 4, 5] + integer(int64) :: i641(5) = [1, 2, 3, 4, 5] + + integer(int32), allocatable :: i32(:,:), i323(:, :, :) + integer(int64), allocatable :: i64(:,:), i643(:, :, :) + + real(sp) :: s1(5) = [1.0_sp, 2.0_sp, 3.0_sp, 4.0_sp, 5.0_sp] + real(dp) :: d1(5) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp] + + real(sp), allocatable :: s(:, :), s3(:, :, :) + real(dp), allocatable :: d3(:, :, :) + real(dp) :: d(4, 3)= reshape([1._dp, 3._dp, 5._dp, 7._dp,& + 2._dp, 4._dp, 6._dp, 8._dp,& + 9._dp, 10._dp, 11._dp, 12._dp], [4, 3]) + + !sp + !1dim + print*,' test_sp_1dim' + call assert( abs(var(s1) - 2.5) < sptol) + call assert( abs(var(s1, dim=1) - 2.5) < sptol) + + print*,' test_sp_1dim_mask' + call assert( isnan(var(s1, .false.))) + call assert( isnan(var(s1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(s1, s1 < 5) - 5./3.) < sptol) + call assert( abs(var(s1, 1, s1 < 5) - 5./3.) < sptol) + + + !2dim + print*,' test_sp_2dim' + s = d + call assert( abs(var(s) - 13) < sptol) + call assert( all( abs( var(s, 1) - [20. / 3., 20. / 3., 5. / 3.]) < sptol)) + call assert( all( abs( var(s, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]) < sptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(s, .false.))) + call assert( any(isnan(var(s, 1, .false.)))) + call assert( any(isnan(var(s, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(s, s < 11) - 27.5 / 3.) < sptol) + call assert( all( abs( var(s, 1, s < 11) - [20. / 3., 20. / 3., 0.5]) < sptol)) + call assert( all( abs( var(s, 2, s < 11) - [19.0, 43. / 3., 0.5 , 0.5]) < sptol)) + + + !3dim + allocate(s3(size(s,1),size(s,2),3)) + s3(:,:,1)=s; + s3(:,:,2)=s*2; + s3(:,:,3)=s*4; + + print*,' test_sp_3dim' + call assert( abs(var(s3) - 153.4) < sptol) + call assert( all( abs( var(s3, 1) -& + reshape([20. / 3., 20. / 3., 5. / 3.,& + 4* 20. / 3., 4* 20. / 3., 4* 5. / 3.,& + 16* 20. / 3., 16* 20. / 3., 16* 5. / 3.],& + [size(s3,2), size(s3,3)]))& + < sptol)) + call assert( all( abs( var(s3, 2) -& + reshape([19.0, 43. / 3., 31. / 3. , 7.0,& + 4* 19.0, 4* 43. / 3., 4* 31. / 3. , 4* 7.0,& + 16* 19.0, 16* 43. / 3., 16* 31. / 3. , 16* 7.0],& + [size(s3,1), size(s3,3)] ))& + < sptol)) + call assert( all(abs( var(s3, 3) -& + reshape([ 7./3., 21., 175./3.,& + 343./3., 28./3., 112./3.,& + 84., 448./3., 189.,& + 700./3., 847./3., 336.], [size(s3,1), size(s3,2)] ))& + < sptol)) + + print*,' test_sp_3dim_mask' + call assert( isnan(var(s3, .false.))) + call assert( any(isnan(var(s3, 1, .false.)))) + call assert( any(isnan(var(s3, 2, .false.)))) + call assert( any(isnan(var(s3, 3, .false.)))) + + print*,' test_sp_3dim_mask_array' + call assert( abs(var(s3, s3 < 11) - 8.2205877_sp) < sptol) + call assert( all( abs( var(s3, 1, s3 < 25) -& + reshape([20./3., 20./3., 5./3., 80./3., 80./3., 20./3., 64., 64., 0.],& + [size(s3, 2), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 2, s3 < 25) -& + reshape([19., 43./3., 31./3., 7.,& + 4*19., 4*43./3., 4*31./3., 4*7.,& + 8., 8., 8., 0.],& + [size(s3, 1), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 3, s3 < 25) -& + reshape([ 7./3., 21., 175./3.,& + 24.5, 28./3., 112./3.,& + 84., 32., 40.5,& + 50., 60.5, 72.], [size(s3,1), size(s3,2)] ))& + < sptol )) + + + !dp + !1dim + print*,' test_dp_1dim' + call assert( abs(var(d1) - 2.5) < dptol) + call assert( abs(var(d1, 1) - 2.5) < dptol) + + print*,' test_dp_1dim_mask' + call assert( isnan(var(d1, .false.))) + call assert( isnan(var(d1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_dp_2dim' + call assert( abs(var(d) - 13) < dptol) + call assert( all( abs( var(d,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(d,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_dp_2dim_mask' + call assert( isnan(var(d, .false.))) + call assert( any(isnan(var(d, 1, .false.)))) + call assert( any(isnan(var(d, 2, .false.)))) + + print*,' test_dp_2dim_mask_array' + call assert( abs(var(d, d < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(d, 1, d < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(d, 2, d < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5 _dp, 0.5_dp]) < dptol)) + + !3dim + allocate(d3(size(d,1),size(d,2),3)) + d3(:,:,1)=d; + d3(:,:,2)=d*2; + d3(:,:,3)=d*4; + + print*,' test_dp_3dim' + call assert( abs(var(d3) - 153.4_dp) < dptol) + call assert( all( abs( var(d3, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(d3,2), size(d3,3)]))& + < dptol)) + print*,' test_dp_3dim' + call assert( all( abs( var(d3, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(d3,1), size(d3,3)] ))& + < dptol)) + print*,' test_dp_3dim' + call assert( all(abs( var(d3, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(d3,1), size(d3,2)] ))& + < dptol)) + + print*,' test_dp_3dim_mask' + call assert( isnan(var(d3, .false.))) + call assert( any(isnan(var(d3, 1, .false.)))) + call assert( any(isnan(var(d3, 2, .false.)))) + call assert( any(isnan(var(d3, 3, .false.)))) + + print*,' test_dp_3dim_mask_array' + call assert( abs(var(d3, d3 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(d3, 1, d3 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(d3, 2), size(d3, 3)]))& + < dptol )) + call assert( all( abs( var(d3, 2, d3 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(d3, 1), size(d3, 3)]))& + < dptol )) + call assert( all( abs( var(d3, 3, d3 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(d3,1), size(d3,2)] ))& + < dptol )) + + + + !int32 + !1dim + print*,' test_int32_1dim' + call assert( abs(var(i321) - 2.5) < dptol) + call assert( abs(var(i321, 1) - 2.5) < dptol) + + print*,' test_int32_1dim_mask' + call assert( isnan(var(i321, .false.))) + call assert( isnan(var(i321, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_int32_2dim' + i32 = d + call assert( abs(var(i32) - 13) < dptol) + call assert( all( abs( var(i32,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(i32,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_int32_2dim_mask' + call assert( isnan(var(i32, .false.))) + call assert( any(isnan(var(i32, 1, .false.)))) + call assert( any(isnan(var(i32, 2, .false.)))) + + print*,' test_int32_2dim_mask_array' + call assert( abs(var(i32, i32 < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(i32, 1, i32 < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(i32, 2, i32 < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol)) + + !3dim + allocate(i323(size(d,1),size(d,2),3)) + i323(:,:,1)=d; + i323(:,:,2)=d*2; + i323(:,:,3)=d*4; + + print*,' test_int32_3dim' + call assert( abs(var(i323) - 153.4_dp) < dptol) + call assert( all( abs( var(i323, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(i323,2), size(i323,3)]))& + < dptol)) + call assert( all( abs( var(i323, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(i323,1), size(i323,3)] ))& + < dptol)) + call assert( all(abs( var(i323, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(i323,1), size(i323,2)] ))& + < dptol)) + + print*,' test_int32_3dim_mask' + call assert( isnan(var(i323, .false.))) + call assert( any(isnan(var(i323, 1, .false.)))) + call assert( any(isnan(var(i323, 2, .false.)))) + call assert( any(isnan(var(i323, 3, .false.)))) + + print*,' test_int32_3dim_mask_array' + call assert( abs(var(i323, i323 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(i323, 1, i323 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(i323, 2), size(i323, 3)]))& + < dptol )) + call assert( all( abs( var(i323, 2, i323 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(i323, 1), size(i323, 3)]))& + < dptol )) + call assert( all( abs( var(i323, 3, i323 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(i323,1), size(i323,2)] ))& + < dptol )) + + + !int64 + !1dim + print*,' test_int64_1dim' + call assert( abs(var(i641) - 2.5) < dptol) + call assert( abs(var(i641, 1) - 2.5) < dptol) + + print*,' test_int64_1dim_mask' + call assert( isnan(var(i641, .false.))) + call assert( isnan(var(i641, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_int64_2dim' + i64 = d + call assert( abs(var(i64) - 13) < dptol) + call assert( all( abs( var(i64,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(i64,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_int64_2dim_mask' + call assert( isnan(var(i64, .false.))) + call assert( any(isnan(var(i64, 1, .false.)))) + call assert( any(isnan(var(i64, 2, .false.)))) + + print*,' test_int64_2dim_mask_array' + call assert( abs(var(i64, i64 < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(i64, 1, i64 < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(i64, 2, i64 < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5 _dp, 0.5_dp]) < dptol)) + + !3dim + allocate(i643(size(d,1),size(d,2),3)) + i643(:,:,1)=d; + i643(:,:,2)=d*2; + i643(:,:,3)=d*4; + + print*,' test_int32_3dim' + call assert( abs(var(i643) - 153.4_dp) < dptol) + call assert( all( abs( var(i643, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(i643,2), size(i643,3)]))& + < dptol)) + call assert( all( abs( var(i643, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(i643,1), size(i643,3)] ))& + < dptol)) + call assert( all(abs( var(i643, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(i643,1), size(i643,2)] ))& + < dptol)) + + print*,' test_int32_3dim_mask' + call assert( isnan(var(i643, .false.))) + call assert( any(isnan(var(i643, 1, .false.)))) + call assert( any(isnan(var(i643, 2, .false.)))) + call assert( any(isnan(var(i643, 3, .false.)))) + + print*,' test_int32_3dim_mask_array' + call assert( abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(i643, 1, i643 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(i643, 2), size(i643, 3)]))& + < dptol )) + call assert( all( abs( var(i643, 2, i643 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(i643, 1), size(i643, 3)]))& + < dptol )) + call assert( all( abs( var(i643, 3, i643 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(i643,1), size(i643,2)] ))& + < dptol )) + +end program From 044abc55653013202c98025fd1ef34d8a23d2adc Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Thu, 6 Feb 2020 00:25:15 +0100 Subject: [PATCH 2/8] varaince_dev: update var modules --- src/common.fypp | 14 ++++++ src/stdlib_experimental_stats.fypp | 48 ++++++++++++-------- src/stdlib_experimental_stats_var.fypp | 63 +++++++++++++++----------- 3 files changed, 79 insertions(+), 46 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index b3e45d3d4..f9ade3112 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -113,6 +113,20 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ $:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix) #:enddef +#! Generates Fortran expressions. +#! +#! Args: +#! varname (str): Name of the variable to be used as origin +#! varname1 (str): Name of the variable to be used instead of varname +#! origrank (int): Rank of the original variable +#! dim (int): Index of the used expression varname1 +#! +#! Returns: +#! Shape expression enclosed in braces, except for the index dim. +#! +#! E.g., (:, :, :, i, :, :) +#! + #:def rankindice(varname, varname1, origrank, dim) #:assert origrank > 0 #:if origrank > 0 diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index ab0582a4d..1958ebe92 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -107,91 +107,99 @@ module stdlib_experimental_stats interface var - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_all_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_all_${k1}$_dp(x, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_all_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res) + #:set RName = rname("var",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 var_${rank}$_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_${k1}$_dp(x, dim, mask) result(res) + #:set RName = rname("var",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 var_${rank}$_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_mask_all_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_mask_all_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_mask_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res) + #:set RName = rname("var_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 var_${rank}$_mask_${k1}$_dp + end function ${RName}$ #:endfor #:endfor diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 675767a0a..b5b37508a 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -1,8 +1,6 @@ #:include "common.fypp" - #:set RANKS = range(1, MAXRANK + 1) - - +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES submodule (stdlib_experimental_stats) stdlib_experimental_stats_var use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan @@ -12,9 +10,10 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_var contains - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res) + #:set RName = rname("var_all",rank, t1, k1) + module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask ${t1}$ :: res @@ -22,7 +21,7 @@ contains ${t1}$ :: n, mean if (.not.optval(mask, .true.)) then - res = ieee_value(res, ieee_quiet_nan) + res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan) return end if @@ -31,14 +30,15 @@ contains res = sum((x - mean)**2) / (n - 1._${k1}$) - end function var_${rank}$_all_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_all_${k1}$_dp(x, mask) result(res) + #:set RName = rname("var_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 @@ -55,14 +55,15 @@ contains res = sum((real(x, dp) - mean)**2) / (n - 1._dp) - end function var_${rank}$_all_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res) + #:set RName = rname("var",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 @@ -73,7 +74,7 @@ contains ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ if (.not.optval(mask, .true.)) then - res = ieee_value(res, ieee_quiet_nan) + res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan) return end if @@ -92,14 +93,15 @@ contains end select res = res / (n - 1._${k1}$) - end function var_${rank}$_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_${k1}$_dp(x, dim, mask) result(res) + #:set RName = rname("var",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 @@ -129,14 +131,15 @@ contains end select res = res / (n - 1._dp) - end function var_${rank}$_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res) + #:set RName = rname("var_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 @@ -148,14 +151,15 @@ contains res = sum((x - mean)**2, mask) / (n - 1._${k1}$) - end function var_${rank}$_mask_all_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res) + #:set RName = rname("var_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 @@ -167,14 +171,15 @@ contains res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp) - end function var_${rank}$_mask_all_${k1}$_dp + end function ${RName}$ #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res) + #:set RName = rname("var_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)}$ @@ -192,7 +197,12 @@ contains mean = sum(x, dim, mask) / n do i = 1, size(x, dim) res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& - 0._${k1}$, mask${rankindice(':', 'i', rank, fi)}$) + #:if t1[0] == 'r' + 0._${k1}$,& + #:else + cmplx(0._${k1}$, 0._${k1}$, ${k1}$),& + #:endif + mask${rankindice(':', 'i', rank, fi)}$) end do #:endfor case default @@ -200,14 +210,15 @@ contains end select res = res / (n - 1._${k1}$) - end function var_${rank}$_mask_${k1}$_${k1}$ + end function ${RName}$ #:endfor #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS - module function var_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res) + #:set RName = rname("var_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)}$ @@ -233,7 +244,7 @@ contains end select res = res / (n - 1._dp) - end function var_${rank}$_mask_${k1}$_dp + end function ${RName}$ #:endfor #:endfor From d77b6e998073ba020e61cf0134078b84db5accc8 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Thu, 6 Feb 2020 00:43:41 +0100 Subject: [PATCH 3/8] variance_dev: update spec var --- src/stdlib_experimental_stats.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index 00f96d6df..f8fd1fb09 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -69,7 +69,7 @@ The variance is defined as the best unbiased estimator and is computed as: ### Arguments -`array`: Shall be an array of type `integer`, or `real`. +`array`: Shall be an array of type `integer`, `real`, or `complex`. `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`. @@ -77,7 +77,7 @@ The variance is defined as the best unbiased estimator and is computed as: ### Return value -If `array` is of type `real`, the result is of the same type as `array`. +If `array` is of type `real` or `complex`, the result is of the same type as `array`. If `array` is of type `integer`, the result is of type `double precision`. If `dim` is absent, a scalar with the variance 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 `ar From baabfc8a5e4133fb26616c7caffbae86e36e629a Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Fri, 7 Feb 2020 21:51:36 +0100 Subject: [PATCH 4/8] variance_dev: changed ieee_value() as proposed --- src/stdlib_experimental_stats_var.fypp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index b5b37508a..0153618b0 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -21,7 +21,7 @@ contains ${t1}$ :: n, mean if (.not.optval(mask, .true.)) then - res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan) + res = ieee_value(1._${k1}$, ieee_quiet_nan) return end if @@ -46,7 +46,7 @@ contains real(dp) :: n, mean if (.not.optval(mask, .true.)) then - res = ieee_value(res, ieee_quiet_nan) + res = ieee_value(1._dp, ieee_quiet_nan) return end if @@ -74,7 +74,7 @@ contains ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ if (.not.optval(mask, .true.)) then - res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan) + res = ieee_value(1._${k1}$, ieee_quiet_nan) return end if @@ -112,7 +112,7 @@ contains real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ if (.not.optval(mask, .true.)) then - res = ieee_value(res, ieee_quiet_nan) + res = ieee_value(1._dp, ieee_quiet_nan) return end if From 9b191545e2da55bb0efa371f5c69435b652752e1 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Fri, 7 Feb 2020 22:01:18 +0100 Subject: [PATCH 5/8] variance_dev: remove support of complex because it was wrong --- src/stdlib_experimental_stats_var.fypp | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 0153618b0..2db92a247 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -1,6 +1,5 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) -#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES submodule (stdlib_experimental_stats) stdlib_experimental_stats_var use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan @@ -10,7 +9,7 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_var contains - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1) module function ${RName}$(x, mask) result(res) @@ -60,7 +59,7 @@ contains #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",rank, t1, k1) module function ${RName}$(x, dim, mask) result(res) @@ -136,7 +135,7 @@ contains #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask_all",rank, t1, k1) module function ${RName}$(x, mask) result(res) @@ -176,7 +175,7 @@ contains #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask",rank, t1, k1) module function ${RName}$(x, dim, mask) result(res) @@ -197,11 +196,7 @@ contains mean = sum(x, dim, mask) / n do i = 1, size(x, dim) res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& - #:if t1[0] == 'r' 0._${k1}$,& - #:else - cmplx(0._${k1}$, 0._${k1}$, ${k1}$),& - #:endif mask${rankindice(':', 'i', rank, fi)}$) end do #:endfor From da90a891d38176ff8e1c095df5eaf0a1935cd4e4 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Fri, 7 Feb 2020 23:10:52 +0100 Subject: [PATCH 6/8] variance_dev: addition of variance of complex as (var(real(x)) + var(aimag(x)) --- src/stdlib_experimental_stats.fypp | 54 ++++++++-- src/stdlib_experimental_stats_var.fypp | 142 +++++++++++++++++++++++++ src/tests/stats/test_var.f90 | 89 +++++++++++++++- 3 files changed, 278 insertions(+), 7 deletions(-) diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index 1958ebe92..b20bb5b45 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -107,7 +107,7 @@ module stdlib_experimental_stats interface var - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1) module function ${RName}$(x, mask) result(res) @@ -118,6 +118,17 @@ module stdlib_experimental_stats #:endfor #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(${k1}$) :: res + end function ${RName}$ + #:endfor + #:endfor + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1, 'dp') @@ -129,7 +140,7 @@ module stdlib_experimental_stats #:endfor #:endfor - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",rank, t1, k1) module function ${RName}$(x, dim, mask) result(res) @@ -141,6 +152,18 @@ module stdlib_experimental_stats #:endfor #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(${k1}$) :: 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("var",rank, t1, k1, 'dp') @@ -153,8 +176,7 @@ module stdlib_experimental_stats #:endfor #:endfor - - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask_all",rank, t1, k1) module function ${RName}$(x, mask) result(res) @@ -165,6 +187,16 @@ module stdlib_experimental_stats #:endfor #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res + end function ${RName}$ + #:endfor + #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS @@ -177,8 +209,7 @@ module stdlib_experimental_stats #:endfor #:endfor - - #:for k1, t1 in RC_KINDS_TYPES + #:for k1, t1 in REAL_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask",rank, t1, k1) module function ${RName}$(x, dim, mask) result(res) @@ -190,6 +221,17 @@ module stdlib_experimental_stats #:endfor #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1, 'r'+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)}$ + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 2db92a247..4056c3667 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -34,6 +34,36 @@ contains #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(${k1}$) :: res + + real(${k1}$) :: n, mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + n = real(size(x, kind = int64), ${k1}$) + + !real part + mean = sum(real(x)) / n + res = sum((real(x) - mean)**2) + + !imaginary part + mean = sum(aimag(x)) / n + res = (res + sum((aimag(x) - mean)**2)) / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1, 'dp') @@ -97,6 +127,50 @@ contains #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(${k1}$) :: n + real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(size(x, dim), ${k1}$) + !real part + mean = sum(real(x), dim) / n + do i = 1, size(x, dim) + res = res + (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2 + end do + !imaginary part + mean = sum(aimag(x), dim) / n + do i = 1, size(x, dim) + res = res + (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2 + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",rank, t1, k1, 'dp') @@ -155,6 +229,31 @@ contains #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res + + real(${k1}$) :: n, mean + + n = real(count(mask, kind = int64), ${k1}$) + + !real part + mean = sum(real(x), mask) / n + res = sum((real(x) - mean)**2, mask) + + !imaginary part + mean = sum(aimag(x), mask) / n + res = (res + sum((aimag(x) - mean)**2, mask)) / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask_all",rank, t1, k1, 'dp') @@ -210,6 +309,49 @@ contains #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1, 'r'+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)}$ + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$ + real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(count(mask, dim), ${k1}$) + !real part + mean = sum(real(x), dim, mask) / n + do i = 1, size(x, dim) + res = res + merge( (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,& + 0._${k1}$,& + mask${rankindice(':', 'i', rank, fi)}$) + end do + !imaginary part + mean = sum(aimag(x), dim, mask) / n + do i = 1, size(x, dim) + res = res + merge( (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,& + 0._${k1}$,& + mask${rankindice(':', 'i', rank, fi)}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask",rank, t1, k1, 'dp') diff --git a/src/tests/stats/test_var.f90 b/src/tests/stats/test_var.f90 index 32481615d..868c7a194 100644 --- a/src/tests/stats/test_var.f90 +++ b/src/tests/stats/test_var.f90 @@ -20,10 +20,25 @@ program test_var real(sp), allocatable :: s(:, :), s3(:, :, :) real(dp), allocatable :: d3(:, :, :) - real(dp) :: d(4, 3)= reshape([1._dp, 3._dp, 5._dp, 7._dp,& + real(dp) :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 7._dp,& 2._dp, 4._dp, 6._dp, 8._dp,& 9._dp, 10._dp, 11._dp, 12._dp], [4, 3]) + + complex(sp) :: cs1(5) = [ cmplx(0.57706_sp, 0.00000_sp),& + cmplx(0.00000_sp, 1.44065_sp),& + cmplx(1.26401_sp, 0.00000_sp),& + cmplx(0.00000_sp, 0.88833_sp),& + cmplx(1.14352_sp, 0.00000_sp)] + complex(dp) :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp),& + cmplx(0.00000_dp, 1.44065_dp),& + cmplx(1.26401_dp, 0.00000_dp),& + cmplx(0.00000_dp, 0.88833_dp),& + cmplx(1.14352_dp, 0.00000_dp)] + complex(sp) :: cs(5,3) + complex(dp) :: cd(5,3) + + !sp !1dim print*,' test_sp_1dim' @@ -388,4 +403,76 @@ program test_var [size(i643,1), size(i643,2)] ))& < dptol )) + !csp + !1dim + print*,' test_csp_1dim' + call assert( abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol) + call assert( abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol) + + print*,' test_csp_1dim_mask' + call assert( isnan(var(cs1, .false.))) + call assert( isnan(var(cs1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol) + call assert( abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol) + + !2dim + cs(:,1) = cs1 + cs(:,2) = cs1*3_sp + cs(:,3) = cs1*1.5_sp + + print*,' test_sp_2dim' + print*,var(cs),(var(real(cs)) + var(aimag(cs))) + print*,var(cs)-(var(real(cs)) + var(aimag(cs))),sptol + call assert( abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol) + call assert( all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol)) + call assert( all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(cs, .false.))) + call assert( any(isnan(var(cs, 1, .false.)))) + call assert( any(isnan(var(cs, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol) + call assert( all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol)) + call assert( all( abs( var(cs, 2, aimag(cs) == 0) - var(real(cs), 2, aimag(cs) == 0)) < sptol)) + + !cdp + !1dim + print*,' test_cdp_1dim' + call assert( abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol) + call assert( abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol) + + print*,' test_cdp_1dim_mask' + call assert( isnan(var(cd1, .false.))) + call assert( isnan(var(cd1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol) + call assert( abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol) + + !2dim + cd(:,1) = cd1 + cd(:,2) = cd1*3_sp + cd(:,3) = cd1*1.5_sp + + print*,' test_sp_2dim' + print*,var(cd),(var(real(cd)) + var(aimag(cd))) + print*,var(cd)-(var(real(cd)) + var(aimag(cd))),dptol + call assert( abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol) + call assert( all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol)) + call assert( all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(cd, .false.))) + call assert( any(isnan(var(cd, 1, .false.)))) + call assert( any(isnan(var(cd, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol) + call assert( all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol)) + call assert( all( abs( var(cd, 2, aimag(cd) == 0) - var(real(cd), 2, aimag(cd) == 0)) < dptol)) + end program From 2a0182aa42b9a501b76638e64968b06f548e46bc Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Fri, 7 Feb 2020 23:54:10 +0100 Subject: [PATCH 7/8] variance_dev: use fypp to avoid abs in real functions --- src/stdlib_experimental_stats.fypp | 54 +------ src/stdlib_experimental_stats_var.fypp | 193 +++++-------------------- 2 files changed, 39 insertions(+), 208 deletions(-) diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index b20bb5b45..ba5e7ad63 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -107,20 +107,9 @@ module stdlib_experimental_stats interface var - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_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 CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_all",rank, t1, k1, 'r'+k1) module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask @@ -140,21 +129,9 @@ module stdlib_experimental_stats #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",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 CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var",rank, t1, k1, 'r'+k1) module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim @@ -176,20 +153,9 @@ module stdlib_experimental_stats #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_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 CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1) module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ @@ -209,21 +175,9 @@ module stdlib_experimental_stats #:endfor #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_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 CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_mask",rank, t1, k1, 'r'+k1) module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 4056c3667..8f210789e 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -1,5 +1,6 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES submodule (stdlib_experimental_stats) stdlib_experimental_stats_var use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan @@ -9,40 +10,16 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_var contains - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1) - module function ${RName}$(x, mask) result(res) - ${t1}$, intent(in) :: x${ranksuffix(rank)}$ - logical, intent(in), optional :: mask - ${t1}$ :: res - - ${t1}$ :: n, mean - - if (.not.optval(mask, .true.)) then - res = ieee_value(1._${k1}$, ieee_quiet_nan) - return - end if - - n = real(size(x, kind = int64), ${k1}$) - mean = sum(x) / n - - res = sum((x - mean)**2) / (n - 1._${k1}$) - - end function ${RName}$ - #:endfor - #:endfor - - - #:for k1, t1 in CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_all",rank, t1, k1, 'r'+k1) module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask real(${k1}$) :: res - real(${k1}$) :: n, mean + real(${k1}$) :: n + ${t1}$ :: mean if (.not.optval(mask, .true.)) then res = ieee_value(1._${k1}$, ieee_quiet_nan) @@ -50,14 +27,13 @@ contains end if n = real(size(x, kind = int64), ${k1}$) + mean = sum(x) / n - !real part - mean = sum(real(x)) / n - res = sum((real(x) - mean)**2) - - !imaginary part - mean = sum(aimag(x)) / n - res = (res + sum((aimag(x) - mean)**2)) / (n - 1._${k1}$) + #:if t1[0] == 'r' + res = sum((x - mean)**2) / (n - 1._${k1}$) + #:else + res = sum(abs(x - mean)**2) / (n - 1._${k1}$) + #:endif end function ${RName}$ #:endfor @@ -89,47 +65,9 @@ contains #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",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 :: i - ${t1}$ :: n - ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ - - if (.not.optval(mask, .true.)) then - res = ieee_value(1._${k1}$, ieee_quiet_nan) - return - end if - - res = 0._${k1}$ - select case(dim) - #:for fi in range(1, rank+1) - case(${fi}$) - n = real(size(x, dim), ${k1}$) - mean = sum(x, dim) / n - do i = 1, size(x, dim) - res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2 - end do - #:endfor - case default - call error_stop("ERROR (mean): wrong dimension") - end select - res = res / (n - 1._${k1}$) - - end function ${RName}$ - #:endfor - #:endfor - - - #:for k1, t1 in CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var",rank, t1, k1, 'r'+k1) module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim @@ -138,7 +76,7 @@ contains integer :: i real(${k1}$) :: n - real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$ + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ if (.not.optval(mask, .true.)) then res = ieee_value(1._${k1}$, ieee_quiet_nan) @@ -150,15 +88,13 @@ contains #:for fi in range(1, rank+1) case(${fi}$) n = real(size(x, dim), ${k1}$) - !real part - mean = sum(real(x), dim) / n - do i = 1, size(x, dim) - res = res + (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2 - end do - !imaginary part - mean = sum(aimag(x), dim) / n + mean = sum(x, dim) / n do i = 1, size(x, dim) - res = res + (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2 + #:if t1[0] == 'r' + res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2 + #:else + res = res + abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2 + #:endif end do #:endfor case default @@ -209,45 +145,25 @@ contains #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_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 - - ${t1}$ :: n, mean - - n = real(count(mask, kind = int64), ${k1}$) - mean = sum(x, mask) / n - - res = sum((x - mean)**2, mask) / (n - 1._${k1}$) - - end function ${RName}$ - #:endfor - #:endfor - - - #:for k1, t1 in CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1) module function ${RName}$(x, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ real(${k1}$) :: res - real(${k1}$) :: n, mean + real(${k1}$) :: n + ${t1}$ :: mean n = real(count(mask, kind = int64), ${k1}$) + mean = sum(x, mask) / n - !real part - mean = sum(real(x), mask) / n - res = sum((real(x) - mean)**2, mask) - - !imaginary part - mean = sum(aimag(x), mask) / n - res = (res + sum((aimag(x) - mean)**2, mask)) / (n - 1._${k1}$) + #:if t1[0] == 'r' + res = sum((x - mean)**2, mask) / (n - 1._${k1}$) + #:else + res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$) + #:endif end function ${RName}$ #:endfor @@ -274,44 +190,9 @@ contains #:endfor - #:for k1, t1 in REAL_KINDS_TYPES + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_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 :: i - ${t1}$ :: n${reduced_shape('x', rank, 'dim')}$ - ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ - - res = 0._${k1}$ - select case(dim) - #:for fi in range(1, rank+1) - case(${fi}$) - n = real(count(mask, dim), ${k1}$) - mean = sum(x, dim, mask) / n - do i = 1, size(x, dim) - res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& - 0._${k1}$,& - mask${rankindice(':', 'i', rank, fi)}$) - end do - #:endfor - case default - call error_stop("ERROR (mean): wrong dimension") - end select - res = res / (n - 1._${k1}$) - - end function ${RName}$ - #:endfor - #:endfor - - - #:for k1, t1 in CMPLX_KINDS_TYPES - #:for rank in RANKS - #:set RName = rname("var_mask",rank, t1, k1, 'r'+k1) module function ${RName}$(x, dim, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim @@ -320,24 +201,20 @@ contains integer :: i real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$ - real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$ + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ res = 0._${k1}$ select case(dim) #:for fi in range(1, rank+1) case(${fi}$) n = real(count(mask, dim), ${k1}$) - !real part - mean = sum(real(x), dim, mask) / n - do i = 1, size(x, dim) - res = res + merge( (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,& - 0._${k1}$,& - mask${rankindice(':', 'i', rank, fi)}$) - end do - !imaginary part - mean = sum(aimag(x), dim, mask) / n + mean = sum(x, dim, mask) / n do i = 1, size(x, dim) - res = res + merge( (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,& + #:if t1[0] == 'r' + res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& + #:else + res = res + merge( abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& + #:endif 0._${k1}$,& mask${rankindice(':', 'i', rank, fi)}$) end do From 01e897cb67de961949e211e57bc38494caedb2bd Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Mon, 10 Feb 2020 16:47:47 +0100 Subject: [PATCH 8/8] variance_dev:suggestions by @aradi --- src/common.fypp | 42 ++++++++++++-------------- src/stdlib_experimental_stats_var.fypp | 16 +++++----- 2 files changed, 27 insertions(+), 31 deletions(-) diff --git a/src/common.fypp b/src/common.fypp index f9ade3112..0dbd3e7d0 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -97,6 +97,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endif #:enddef + #! Generates a routine name from a generic name, rank, type and kind #! #! Args: @@ -113,33 +114,28 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ $:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix) #:enddef -#! Generates Fortran expressions. + +#! Generates an array rank suffix for subarrays reducing the dimension #! #! Args: -#! varname (str): Name of the variable to be used as origin -#! varname1 (str): Name of the variable to be used instead of varname -#! origrank (int): Rank of the original variable -#! dim (int): Index of the used expression varname1 +#! rank (int): Rank of the original variable +#! selectors (array): Dimension and name of the variable(s) #! #! Returns: -#! Shape expression enclosed in braces, except for the index dim. -#! -#! E.g., (:, :, :, i, :, :) -#! - -#:def rankindice(varname, varname1, origrank, dim) - #:assert origrank > 0 - #:if origrank > 0 - #:call join_lines(joinstr=", ", prefix="(", suffix=")") - #:for i in range(1, origrank+1) - #:if i == dim - ${varname1}$ - #:else - ${varname}$ - #:endif - #:endfor - #:endcall - #:endif +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! select_subarray(5 , [(4, 'i'), (5, 'j')])}$ +#! -> (:, :, :, i, j) +#! +#:def select_subarray(rank, selectors) + #:assert rank > 0 + #:set seldict = dict(selectors) + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + $:seldict.get(i, ":") + #:endfor + #:endcall #:enddef #:endmute diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 8f210789e..eaa5bdcdd 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -91,9 +91,9 @@ contains mean = sum(x, dim) / n do i = 1, size(x, dim) #:if t1[0] == 'r' - res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2 + res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2 #:else - res = res + abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2 + res = res + abs(x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2 #:endif end do #:endfor @@ -132,7 +132,7 @@ contains n = real(size(x, dim), dp) mean = sum(real(x, dp), dim) / n do i = 1, size(x, dim) - res = res + (real(x${rankindice(':', 'i', rank, fi )}$, dp) - mean)**2 + res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2 end do #:endfor case default @@ -211,12 +211,12 @@ contains mean = sum(x, dim, mask) / n do i = 1, size(x, dim) #:if t1[0] == 'r' - res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& + res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,& #:else - res = res + merge( abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2,& + res = res + merge( abs(x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,& #:endif 0._${k1}$,& - mask${rankindice(':', 'i', rank, fi)}$) + mask${select_subarray(rank, [(fi, 'i')])}$) end do #:endfor case default @@ -249,8 +249,8 @@ contains n = real(count(mask, dim), dp) mean = sum(real(x, dp), dim, mask) / n do i = 1, size(x, dim) - res = res + merge((real(x${rankindice(':', 'i', rank, fi )}$, dp) - mean)**2,& - 0._dp, mask${rankindice(':', 'i', rank, fi)}$) + res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2,& + 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) end do #:endfor case default