From 50367ff510e3a871669c3727d9b3f1d14fd212e9 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Sat, 22 Feb 2020 22:10:09 +0100 Subject: [PATCH 1/4] Addition of corrected in API of var --- src/stdlib_experimental_stats.fypp | 24 +++-- src/stdlib_experimental_stats.md | 28 ++++-- src/stdlib_experimental_stats_var.fypp | 53 ++++++---- src/tests/stats/CMakeLists.txt | 1 + src/tests/stats/test_var.f90 | 88 +++++++---------- src/tests/stats/test_varn.f90 | 129 +++++++++++++++++++++++++ 6 files changed, 237 insertions(+), 86 deletions(-) create mode 100644 src/tests/stats/test_varn.f90 diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index ba5e7ad63..c74f6ff6f 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -110,9 +110,10 @@ module stdlib_experimental_stats #: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) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(${k1}$) :: res end function ${RName}$ #:endfor @@ -121,9 +122,10 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1, 'dp') - module function ${RName}$(x, mask) result(res) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(dp) :: res end function ${RName}$ #:endfor @@ -132,10 +134,11 @@ module stdlib_experimental_stats #: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) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor @@ -144,10 +147,11 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(dp) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor @@ -156,9 +160,10 @@ module stdlib_experimental_stats #: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) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(${k1}$) :: res end function ${RName}$ #:endfor @@ -167,9 +172,10 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask_all",rank, t1, k1, 'dp') - module function ${RName}$(x, mask) result(res) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(dp) :: res end function ${RName}$ #:endfor @@ -178,10 +184,11 @@ module stdlib_experimental_stats #: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) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor @@ -190,10 +197,11 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask",rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(dp) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index f8fd1fb09..85b07b58f 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -55,17 +55,22 @@ end program demo_mean 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: +Per default, 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 ``` +where n is the number of elements. + +The use of the term `n-1` for scaling is called Bessel 's correction. The scaling can be changed with the logical argument `corrected`. If `corrected` is `.false.`, then the sum is scaled with `n`, otherwise with `n-1`. + + ### Syntax -`result = var(array [, mask])` +`result = var(array [, mask [, corrected]])` -`result = var(array, dim [, mask])` +`result = var(array, dim [, mask [, corrected]])` ### Arguments @@ -75,16 +80,19 @@ The variance is defined as the best unbiased estimator and is computed as: `mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`. +`corrected` (optional): Shall be a scalar of type `logical`. If `corrected` is `.true.` (default value), the sum is scaled with `n-1`. If `corrected` is `.false.`, then the sum is scaled with `n`. + ### Return value -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 `array` is of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`. +If `array` is of type `integer`, the result is of type `real(dp)`. -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 `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 `array` 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`. +If the variance is computed with only one single element, then the result is IEEE `NaN` if `corrected` is `.true.` and is `0.` if `corrected` is `.false.`. + ### Example ```fortran @@ -93,10 +101,14 @@ program demo_var implicit none real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ] print *, var(x) !returns 3.5 + print *, var(x, corrected = .false.) !returns 2.9167 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] + reshape(x, [ 2, 3 ] ) > 3.) !returns [NaN, NaN, 0.5] + print *, var( reshape(x, [ 2, 3 ] ), 1,& + reshape(x, [ 2, 3 ] ) > 3.,& + corrected=.false.) !returns [NaN, 0., 0.5] end program demo_var ``` diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index eaa5bdcdd..bc66cf810 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -13,9 +13,10 @@ contains #: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) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(${k1}$) :: res real(${k1}$) :: n @@ -30,9 +31,11 @@ contains mean = sum(x) / n #:if t1[0] == 'r' - res = sum((x - mean)**2) / (n - 1._${k1}$) + res = sum((x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& + (optval(corrected, .true.) .and. n > 0._${k1}$))) #:else - res = sum(abs(x - mean)**2) / (n - 1._${k1}$) + res = sum(abs(x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& + (optval(corrected, .true.) .and. n > 0._${k1}$))) #:endif end function ${RName}$ @@ -43,9 +46,10 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_all",rank, t1, k1, 'dp') - module function ${RName}$(x, mask) result(res) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(dp) :: res real(dp) :: n, mean @@ -58,7 +62,8 @@ contains n = real(size(x, kind = int64), dp) mean = sum(real(x, dp)) / n - res = sum((real(x, dp) - mean)**2) / (n - 1._dp) + res = sum((real(x, dp) - mean)**2) / (n - merge(1._dp, 0._dp,& + (optval(corrected, .true.) .and. n > 0._dp))) end function ${RName}$ #:endfor @@ -68,10 +73,11 @@ contains #: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) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ integer :: i @@ -100,7 +106,8 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - 1._${k1}$) + res = res / (n - merge(1._${k1}$, 0._${k1}$,& + (optval(corrected, .true.) .and. n > 0._${k1}$))) end function ${RName}$ #:endfor @@ -110,10 +117,11 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var",rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in), optional :: mask + logical, intent(in), optional :: corrected real(dp) :: res${reduced_shape('x', rank, 'dim')}$ integer :: i @@ -138,7 +146,8 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - 1._dp) + res = res / (n - merge(1._dp, 0._dp,& + (optval(corrected, .true.) .and. n > 0._dp))) end function ${RName}$ #:endfor @@ -148,9 +157,10 @@ contains #: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) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(${k1}$) :: res real(${k1}$) :: n @@ -160,9 +170,11 @@ contains mean = sum(x, mask) / n #:if t1[0] == 'r' - res = sum((x - mean)**2, mask) / (n - 1._${k1}$) + res = sum((x - mean)**2, mask) / (n -& + merge(1._${k1}$, 0._${k1}$, (optval(corrected, .true.)) .and. n > 0._${k1}$)) #:else - res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$) + res = sum(abs(x - mean)**2, mask) / (n - merge(1._${k1}$, 0._${k1}$,& + (optval(corrected, .true.) .and. n > 0._${k1}$))) #:endif end function ${RName}$ @@ -173,9 +185,10 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask_all",rank, t1, k1, 'dp') - module function ${RName}$(x, mask) result(res) + module function ${RName}$(x, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(dp) :: res real(dp) :: n, mean @@ -183,7 +196,8 @@ contains 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) + res = sum((real(x, dp) - mean)**2, mask) / (n - merge(1._dp, 0._dp,& + (optval(corrected, .true.) .and. n > 0._dp))) end function ${RName}$ #:endfor @@ -193,10 +207,11 @@ contains #: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) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ integer :: i @@ -222,7 +237,8 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - 1._${k1}$) + res = res / (n - merge(1._${k1}$, 0._${k1}$,& + (optval(corrected, .true.) .and. n > 0._${k1}$))) end function ${RName}$ #:endfor @@ -232,10 +248,11 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("var_mask",rank, t1, k1, 'dp') - module function ${RName}$(x, dim, mask) result(res) + module function ${RName}$(x, dim, mask, corrected) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: dim logical, intent(in) :: mask${ranksuffix(rank)}$ + logical, intent(in), optional :: corrected real(dp) :: res${reduced_shape('x', rank, 'dim')}$ integer :: i @@ -256,7 +273,7 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - 1._dp) + res = res / (n - merge(1._dp, 0._dp, (optval(corrected, .true.) .and. n > 0._dp))) end function ${RName}$ #:endfor diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 239ac2373..8dc45370c 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,5 +1,6 @@ ADDTEST(mean) ADDTEST(var) +ADDTEST(varn) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/test_var.f90 b/src/tests/stats/test_var.f90 index 868c7a194..fa339d080 100644 --- a/src/tests/stats/test_var.f90 +++ b/src/tests/stats/test_var.f90 @@ -1,8 +1,7 @@ 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 + use stdlib_experimental_stats, only: var implicit none @@ -51,9 +50,10 @@ program test_var print*,' test_sp_1dim_mask_array' call assert( abs(var(s1, s1 < 5) - 5./3.) < sptol) + call assert( isnan((var(s1, s1 < 0.)))) + call assert( isnan((var(s1, s1 == 1.)))) call assert( abs(var(s1, 1, s1 < 5) - 5./3.) < sptol) - !2dim print*,' test_sp_2dim' s = d @@ -107,14 +107,11 @@ program test_var 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.],& + call assert( all( abs( var(s3, 1, s3 < 45) -& + reshape([20./3., 20./3., 5./3., 80./3., 80./3., 20./3.,& + 320./3., 320./3., 16.],& [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( any( isnan( var(s3, 2, s3 < 25)))) call assert( all( abs( var(s3, 3, s3 < 25) -& reshape([ 7./3., 21., 175./3.,& 24.5, 28./3., 112./3.,& @@ -133,8 +130,10 @@ program test_var call assert( isnan(var(d1, .false.))) call assert( isnan(var(d1, 1, .false.))) - print*,' test_sp_1dim_mask_array' + print*,' test_dp_1dim_mask_array' call assert( abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol) + call assert( isnan((var(d1, d1 < 0.)))) + call assert( isnan((var(d1, d1 == 1.)))) call assert( abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol) !2dim @@ -196,18 +195,13 @@ program test_var 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) -& + call assert( all( abs( var(d3, 1, d3 < 45) -& 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],& + 320._dp/3._dp, 320._dp/3._dp, 16._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( any( isnan( var(d3, 2, d3 < 25)))) 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,& @@ -228,8 +222,10 @@ program test_var call assert( isnan(var(i321, .false.))) call assert( isnan(var(i321, 1, .false.))) - print*,' test_sp_1dim_mask_array' + print*,' test_int32_1dim_mask_array' call assert( abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol) + call assert( isnan((var(i321, i321 < 0)))) + call assert( isnan((var(i321, i321 == 1)))) call assert( abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol) !2dim @@ -290,18 +286,13 @@ program test_var 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) -& + call assert( all( abs( var(i323, 1, i323 < 45) -& 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],& + 320._dp/3._dp, 320._dp/3._dp, 16._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( any( isnan( var(i323, 2, i323 < 25)))) 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,& @@ -321,8 +312,10 @@ program test_var call assert( isnan(var(i641, .false.))) call assert( isnan(var(i641, 1, .false.))) - print*,' test_sp_1dim_mask_array' + print*,' test_int641_1dim_mask_array' call assert( abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol) + call assert( isnan((var(i641, i641 < 0)))) + call assert( isnan((var(i641, i641 == 1)))) call assert( abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol) !2dim @@ -381,20 +374,15 @@ program test_var call assert( any(isnan(var(i643, 2, .false.)))) call assert( any(isnan(var(i643, 3, .false.)))) - print*,' test_int32_3dim_mask_array' + print*,' test_int64_3dim_mask_array' call assert( abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol) - call assert( all( abs( var(i643, 1, i643 < 25) -& + call assert( all( abs( var(i643, 1, i643 < 45) -& 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],& + 320._dp/3._dp, 320._dp/3._dp, 16._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( any( isnan( var(i643, 2, i643 < 25)))) 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,& @@ -413,7 +401,7 @@ program test_var call assert( isnan(var(cs1, .false.))) call assert( isnan(var(cs1, 1, .false.))) - print*,' test_sp_1dim_mask_array' + print*,' test_csp_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) @@ -422,22 +410,20 @@ program test_var 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 + print*,' test_csp_2dim' 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' + print*,' test_csp_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' + print*,' test_csp_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)) + call assert( any( isnan( var(cs, 2, aimag(cs) == 0)))) !cdp !1dim @@ -449,7 +435,7 @@ program test_var call assert( isnan(var(cd1, .false.))) call assert( isnan(var(cd1, 1, .false.))) - print*,' test_sp_1dim_mask_array' + print*,' test_cdp_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) @@ -458,21 +444,19 @@ program test_var 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 + print*,' test_cdp_2dim' 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' + print*,' test_cdp_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' + print*,' test_cdp_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)) + call assert( any( isnan( var(cd, 2, aimag(cd) == 0)))) end program diff --git a/src/tests/stats/test_varn.f90 b/src/tests/stats/test_varn.f90 new file mode 100644 index 000000000..5db0bed60 --- /dev/null +++ b/src/tests/stats/test_varn.f90 @@ -0,0 +1,129 @@ +program test_varn + use stdlib_experimental_error, only: assert + use stdlib_experimental_kinds, only: sp, dp, int32, int64 + use stdlib_experimental_stats, only: 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]) + + + 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' + call assert( abs(var(s1, corrected=.false.) - 2.5*(4./5.)) < sptol) + call assert( abs(var(s1, dim=1, corrected=.false.) - 2.5*(4./5.)) < sptol) + + print*,' test_sp_1dim_mask' + call assert( isnan(var(s1, .false., corrected=.false.))) + call assert( isnan(var(s1, 1, .false., corrected=.false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(s1, s1 < 5, corrected=.false.) - 5./4.) < sptol) + call assert( isnan((var(s1, s1 < 0., corrected=.false.)))) + call assert( abs(var(s1, s1 == 1., corrected=.false.)) < sptol) + call assert( abs(var(s1, 1, s1 < 5, corrected=.false.) - 5./4.) < sptol) + + !2dim + print*,' test_sp_2dim' + s = d + call assert( abs(var(s, corrected=.false.) - 13.*11./12.) < sptol) + call assert( all( abs( var(s, 1, corrected=.false.) - [20., 20., 5.]/4.) < sptol)) + call assert( all( abs( var(s, 2, corrected=.false.) -& + [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(s, .false., corrected=.false.))) + call assert( any(isnan(var(s, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(s, 2, .false., corrected=.false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(s, s < 11, corrected=.false.) - 2.75*3.) < sptol) + call assert( all( abs( var(s, 1, s < 11, corrected=.false.) -& + [5., 5., 0.25]) < sptol)) + call assert( all( abs( var(s, 2, s < 11, corrected=.false.) -& + [19.0*2./3., 43./9.*2., 0.25 , 0.25]) < 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, corrected=.false.) - 153.4*35./36.) < sptol) + call assert( all( abs( var(s3, 1, corrected=.false.) -& + 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)])*3./4.)& + < sptol)) + call assert( all( abs( var(s3, 2, corrected=.false.) -& + 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)] )*2./3.)& + < sptol)) + call assert( all(abs( var(s3, 3, corrected=.false.) -& + 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)] )*2./3.)& + < sptol)) + + print*,' test_sp_3dim_mask' + call assert( isnan(var(s3, .false., corrected=.false.))) + call assert( any(isnan(var(s3, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(s3, 2, .false., corrected=.false.)))) + call assert( any(isnan(var(s3, 3, .false., corrected=.false.)))) + + print*,' test_sp_3dim_mask_array' + call assert( abs(var(s3, s3 < 11, corrected=.false.) - 7.73702383_sp) < sptol) + call assert( all( abs( var(s3, 1, s3 < 45, corrected=.false.) -& + reshape([5., 5., 1.25, 20., 20., 5., 80., 80., 32./3.],& + [size(s3, 2), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 2, s3 < 45, corrected=.false.) -& + reshape([ 38./3., 86./9., 6.88888931, 14./3., 152./3.,& + 38.2222214, 27.5555573, 18.6666660, 202.666672,& + 152.888885, 110.222229, 4.& + ],& + [size(s3, 1), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 3, s3 < 45, corrected=.false.) -& + reshape([1.555555, 14., 38.888888, 76.222222, 6.2222222,& + 24.888888, 56., 99.5555, 126., 155.555555, 188.22222, 36.& + ], [size(s3,1), size(s3,2)] ))& + < sptol )) + + +end program From fdfd7f0d9a7f2a563bc194ade7080f6d40a1357a Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Sat, 22 Feb 2020 22:35:53 +0100 Subject: [PATCH 2/4] var_corrected_dev:remove a unnecessary condition --- src/stdlib_experimental_stats_var.fypp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index bc66cf810..7bf5ca77a 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -32,10 +32,10 @@ contains #:if t1[0] == 'r' res = sum((x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& - (optval(corrected, .true.) .and. n > 0._${k1}$))) + optval(corrected, .true.))) #:else res = sum(abs(x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& - (optval(corrected, .true.) .and. n > 0._${k1}$))) + optval(corrected, .true.))) #:endif end function ${RName}$ @@ -63,7 +63,7 @@ contains mean = sum(real(x, dp)) / n res = sum((real(x, dp) - mean)**2) / (n - merge(1._dp, 0._dp,& - (optval(corrected, .true.) .and. n > 0._dp))) + optval(corrected, .true.))) end function ${RName}$ #:endfor @@ -107,7 +107,7 @@ contains call error_stop("ERROR (mean): wrong dimension") end select res = res / (n - merge(1._${k1}$, 0._${k1}$,& - (optval(corrected, .true.) .and. n > 0._${k1}$))) + optval(corrected, .true.))) end function ${RName}$ #:endfor @@ -147,7 +147,7 @@ contains call error_stop("ERROR (mean): wrong dimension") end select res = res / (n - merge(1._dp, 0._dp,& - (optval(corrected, .true.) .and. n > 0._dp))) + optval(corrected, .true.))) end function ${RName}$ #:endfor From 21f68ddbf04762ae1ba4288defb6e5b4c9474ceb Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Sun, 23 Feb 2020 19:55:52 +0100 Subject: [PATCH 3/4] var_corrected_dev: review following comments in #151 --- src/stdlib_experimental_stats_var.fypp | 45 +++++++++++--------------- 1 file changed, 19 insertions(+), 26 deletions(-) diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp index 7bf5ca77a..d80de5c7d 100644 --- a/src/stdlib_experimental_stats_var.fypp +++ b/src/stdlib_experimental_stats_var.fypp @@ -27,15 +27,13 @@ contains return end if - n = real(size(x, kind = int64), ${k1}$) + n = size(x, kind = int64) mean = sum(x) / n #:if t1[0] == 'r' - res = sum((x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& - optval(corrected, .true.))) + res = sum((x - mean)**2) / (n - merge(1, 0 , optval(corrected, .true.))) #:else - res = sum(abs(x - mean)**2) / (n - merge(1._${k1}$, 0._${k1}$,& - optval(corrected, .true.))) + res = sum(abs(x - mean)**2) / (n - merge(1, 0, optval(corrected, .true.))) #:endif end function ${RName}$ @@ -59,11 +57,10 @@ contains return end if - n = real(size(x, kind = int64), dp) + n = size(x, kind = int64) mean = sum(real(x, dp)) / n - res = sum((real(x, dp) - mean)**2) / (n - merge(1._dp, 0._dp,& - optval(corrected, .true.))) + res = sum((real(x, dp) - mean)**2) / (n - merge(1, 0, optval(corrected, .true.))) end function ${RName}$ #:endfor @@ -93,7 +90,7 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - n = real(size(x, dim), ${k1}$) + n = size(x, dim) mean = sum(x, dim) / n do i = 1, size(x, dim) #:if t1[0] == 'r' @@ -106,8 +103,7 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - merge(1._${k1}$, 0._${k1}$,& - optval(corrected, .true.))) + res = res / (n - merge(1, 0, optval(corrected, .true.))) end function ${RName}$ #:endfor @@ -137,7 +133,7 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - n = real(size(x, dim), dp) + n = size(x, dim) mean = sum(real(x, dp), dim) / n do i = 1, size(x, dim) res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2 @@ -146,8 +142,7 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - merge(1._dp, 0._dp,& - optval(corrected, .true.))) + res = res / (n - merge(1, 0, optval(corrected, .true.))) end function ${RName}$ #:endfor @@ -166,16 +161,15 @@ contains real(${k1}$) :: n ${t1}$ :: mean - n = real(count(mask, kind = int64), ${k1}$) + n = count(mask, kind = int64) mean = sum(x, mask) / n #:if t1[0] == 'r' res = sum((x - mean)**2, mask) / (n -& - merge(1._${k1}$, 0._${k1}$, (optval(corrected, .true.)) .and. n > 0._${k1}$)) #:else - res = sum(abs(x - mean)**2, mask) / (n - merge(1._${k1}$, 0._${k1}$,& - (optval(corrected, .true.) .and. n > 0._${k1}$))) + res = sum(abs(x - mean)**2, mask) / (n -& #:endif + merge(1, 0, (optval(corrected, .true.) .and. n > 0))) end function ${RName}$ #:endfor @@ -193,11 +187,11 @@ contains real(dp) :: n, mean - n = real(count(mask, kind = int64), dp) + n = count(mask, kind = int64) mean = sum(real(x, dp), mask) / n - res = sum((real(x, dp) - mean)**2, mask) / (n - merge(1._dp, 0._dp,& - (optval(corrected, .true.) .and. n > 0._dp))) + res = sum((real(x, dp) - mean)**2, mask) / (n -& + merge(1, 0, (optval(corrected, .true.) .and. n > 0))) end function ${RName}$ #:endfor @@ -222,7 +216,7 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - n = real(count(mask, dim), ${k1}$) + n = count(mask, dim) mean = sum(x, dim, mask) / n do i = 1, size(x, dim) #:if t1[0] == 'r' @@ -237,8 +231,7 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - merge(1._${k1}$, 0._${k1}$,& - (optval(corrected, .true.) .and. n > 0._${k1}$))) + res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0))) end function ${RName}$ #:endfor @@ -263,7 +256,7 @@ contains select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - n = real(count(mask, dim), dp) + n = count(mask, dim) mean = sum(real(x, dp), dim, mask) / n do i = 1, size(x, dim) res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2,& @@ -273,7 +266,7 @@ contains case default call error_stop("ERROR (mean): wrong dimension") end select - res = res / (n - merge(1._dp, 0._dp, (optval(corrected, .true.) .and. n > 0._dp))) + res = res / (n - merge(1, 0, (optval(corrected, .true.) .and. n > 0))) end function ${RName}$ #:endfor From 9af45ef098b2f72a527359b16864a0cdcddf9783 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Sun, 23 Feb 2020 21:39:32 +0100 Subject: [PATCH 4/4] var_corrected_dev: addition of tests for dp, int, and complex(dp) for var(..., corrected=.false.) --- src/tests/stats/test_varn.f90 | 245 ++++++++++++++++++++++++++++++++-- 1 file changed, 237 insertions(+), 8 deletions(-) diff --git a/src/tests/stats/test_varn.f90 b/src/tests/stats/test_varn.f90 index 5db0bed60..ea8766bae 100644 --- a/src/tests/stats/test_varn.f90 +++ b/src/tests/stats/test_varn.f90 @@ -9,10 +9,8 @@ program test_varn 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] @@ -24,17 +22,11 @@ program test_varn 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) @@ -125,5 +117,242 @@ program test_varn ], [size(s3,1), size(s3,2)] ))& < sptol )) + !dp + !1dim + print*,' test_dp_1dim' + call assert( abs(var(d1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol) + call assert( abs(var(d1, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol) + + print*,' test_dp_1dim_mask' + call assert( isnan(var(d1, .false., corrected=.false.))) + call assert( isnan(var(d1, 1, .false., corrected=.false.))) + + print*,' test_dp_1dim_mask_array' + call assert( abs(var(d1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol) + call assert( isnan((var(d1, d1 < 0, corrected=.false.)))) + call assert( abs(var(d1, d1 == 1, corrected=.false.)) < dptol) + call assert( abs(var(d1, 1, d1 < 5, corrected=.false.) - 5._dp/4.) < dptol) + + !2dim + print*,' test_dp_2dim' + call assert( abs(var(d, corrected=.false.) - 13._dp*11./12.) < dptol) + call assert( all( abs( var(d, 1, corrected=.false.) - [20., 20., 5.]/4._dp) < dptol)) + call assert( all( abs( var(d, 2, corrected=.false.) -& + [38._dp, 86._dp / 3._dp, 62._dp / 3._dp , 14._dp]/3._dp) < dptol)) + + print*,' test_dp_2dim_mask' + call assert( isnan(var(d, .false., corrected=.false.))) + call assert( any(isnan(var(d, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(d, 2, .false., corrected=.false.)))) + + print*,' test_dp_2dim_mask_array' + call assert( abs(var(d, d < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol) + call assert( all( abs( var(d, 1, d < 11, corrected=.false.) -& + [5._dp, 5._dp, 0.25_dp]) < dptol)) + call assert( all( abs( var(d, 2, d < 11, corrected=.false.) -& + [38._dp/3., 86._dp/9., 0.25_dp , 0.25_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, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol) + call assert( all( abs( var(d3, 1, corrected=.false.) -& + reshape([20._dp , 20._dp , 5._dp ,& + 4* 20._dp , 4* 20._dp , 4* 5._dp ,& + 16* 20._dp , 16* 20._dp , 16* 5._dp ],& + [size(d3,2), size(d3,3)])/4._dp)& + < dptol)) + call assert( all( abs( var(d3, 2, corrected=.false.) -& + reshape([38._dp, 86. / 3._dp, 62. / 3._dp , 14._dp,& + 8* 19._dp, 8* 43. / 3._dp, 8* 31. / 3._dp, 8* 7._dp,& + 32* 19._dp, 32* 43. / 3._dp, 32* 31. / 3._dp, 32* 7._dp],& + [size(d3,1), size(d3,3)] )/3._dp)& + < dptol)) + print*,' test_dp_3dim' + call assert( all(abs( var(d3, 3, corrected=.false.) -& + reshape([ 7._dp/3., 21._dp, 175._dp/3.,& + 343._dp/3., 28._dp/3., 112._dp/3.,& + 84._dp, 448._dp/3., 189._dp,& + 700._dp/3., 847._dp/3., 336._dp],& + [size(d3,1), size(d3,2)] )*2._dp/3._dp)& + < dptol)) + + print*,' test_dp_3dim_mask' + call assert( isnan(var(d3, .false., corrected=.false.))) + call assert( any(isnan(var(d3, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(d3, 2, .false., corrected=.false.)))) + call assert( any(isnan(var(d3, 3, .false., corrected=.false.)))) + + print*,' test_dp_3dim_mask_array' + call assert( abs(var(d3, d3 < 11, corrected=.false.) -& + 7.7370242214532876_dp) < dptol) + call assert( all( abs( var(d3, 1, d3 < 45, corrected=.false.) -& + reshape([5._dp, 5._dp, 1.25_dp, 20._dp, 20._dp, 5._dp,& + 80._dp, 80._dp, 32._dp/3.],& + [size(d3, 2), size(d3, 3)])) < dptol )) + call assert( all( abs( var(d3, 2, d3 < 45, corrected=.false.) -& + reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,& + 344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,& + 1376._dp/9., 992._dp/9., 4._dp& + ],& + [size(d3, 1), size(d3, 3)])) < dptol )) + print*,' test_dp_3dim_mask_array' + call assert( all( abs( var(d3, 3, d3 < 45, corrected=.false.) -& + reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,& + 224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,& + 1694._dp/9., 36._dp& + ], [size(d3,1), size(d3,2)] ))& + < dptol )) + + !int32 + !1dim + print*,' test_int32_1dim' + call assert( abs(var(i321, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol) + call assert( abs(var(i321, dim=1, corrected=.false.) - 2.5_dp*(4._dp/5.)) < dptol) + + print*,' test_int32_1dim_mask' + call assert( isnan(var(i321, .false., corrected=.false.))) + call assert( isnan(var(i321, 1, .false., corrected=.false.))) + + print*,' test_int32_1dim_mask_array' + call assert( abs(var(i321, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol) + call assert( isnan((var(i321, i321 < 0, corrected=.false.)))) + call assert( abs(var(i321, i321 == 1, corrected=.false.)) < dptol) + call assert( abs(var(i321, 1, i321 < 5, corrected=.false.) - 5._dp/4.) < dptol) + + !2dim + i32 = d + print*,' test_int32_2dim' + call assert( abs(var(i32, corrected=.false.) - 13._dp*11./12.) < dptol) + call assert( all( abs( var(i32, 1, corrected=.false.) -& + [20., 20., 5.]/4._dp) < dptol)) + call assert( all( abs( var(i32, 2, corrected=.false.) -& + [38._dp, 86._dp / 3._dp, 62._dp / 3._dp , 14._dp]/3._dp) < dptol)) + + print*,' test_int32_2dim_mask' + call assert( isnan(var(i32, .false., corrected=.false.))) + call assert( any(isnan(var(i32, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(i32, 2, .false., corrected=.false.)))) + + print*,' test_int32_2dim_mask_array' + call assert( abs(var(i32, i32 < 11, corrected=.false.) - 2.75_dp*3._dp) < dptol) + call assert( all( abs( var(i32, 1, i32 < 11, corrected=.false.) -& + [5._dp, 5._dp, 0.25_dp]) < dptol)) + call assert( all( abs( var(i32, 2, i32 < 11, corrected=.false.) -& + [38._dp/3., 86._dp/9., 0.25_dp , 0.25_dp]) < dptol)) + + + !3dim + allocate(i323(size(i32,1),size(i32,2),3)) + i323(:,:,1)=i32; + i323(:,:,2)=i32*2; + i323(:,:,3)=i32*4; + + print*,' test_int32_3dim' + call assert( abs(var(i323, corrected=.false.) - 153.4_dp*35._dp/36._dp) < dptol) + call assert( all( abs( var(i323, 1, corrected=.false.) -& + reshape([20._dp , 20._dp , 5._dp ,& + 4* 20._dp , 4* 20._dp , 4* 5._dp ,& + 16* 20._dp , 16* 20._dp , 16* 5._dp ],& + [size(i323,2), size(i323,3)])/4._dp)& + < dptol)) + call assert( all( abs( var(i323, 2, corrected=.false.) -& + reshape([38._dp, 86. / 3._dp, 62. / 3._dp , 14._dp,& + 8* 19._dp, 8* 43. / 3._dp, 8* 31. / 3._dp, 8* 7._dp,& + 32* 19._dp, 32* 43. / 3._dp, 32* 31. / 3._dp, 32* 7._dp],& + [size(i323,1), size(i323,3)] )/3._dp)& + < dptol)) + print*,' test_int32_3dim' + call assert( all(abs( var(i323, 3, corrected=.false.) -& + reshape([ 7._dp/3., 21._dp, 175._dp/3.,& + 343._dp/3., 28._dp/3., 112._dp/3.,& + 84._dp, 448._dp/3., 189._dp,& + 700._dp/3., 847._dp/3., 336._dp],& + [size(i323,1), size(i323,2)] )*2._dp/3._dp)& + < dptol)) + + print*,' test_int32_3dim_mask' + call assert( isnan(var(i323, .false., corrected=.false.))) + call assert( any(isnan(var(i323, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(i323, 2, .false., corrected=.false.)))) + call assert( any(isnan(var(i323, 3, .false., corrected=.false.)))) + + print*,' test_int32_3dim_mask_array' + call assert( abs(var(i323, i323 < 11, corrected=.false.) -& + 7.7370242214532876_dp) < dptol) + call assert( all( abs( var(i323, 1, i323 < 45, corrected=.false.) -& + reshape([5._dp, 5._dp, 1.25_dp, 20._dp, 20._dp, 5._dp,& + 80._dp, 80._dp, 32._dp/3.],& + [size(i323, 2), size(i323, 3)])) < dptol )) + call assert( all( abs( var(i323, 2, i323 < 45, corrected=.false.) -& + reshape([ 38._dp/3., 86._dp/9., 62._dp/9., 14._dp/3., 152._dp/3.,& + 344._dp/9., 248._dp/9., 168._dp/9., 1824._dp/9.,& + 1376._dp/9., 992._dp/9., 4._dp& + ],& + [size(i323, 1), size(i323, 3)])) < dptol )) + print*,' test_int32_3dim_mask_array' + call assert( all( abs( var(i323, 3, i323 < 45, corrected=.false.) -& + reshape([14._dp/9., 14._dp, 350._dp/9., 686._dp/9., 56._dp/9.,& + 224._dp/9., 56._dp, 896._dp/9., 126._dp, 1400._dp/9.,& + 1694._dp/9., 36._dp& + ], [size(i323,1), size(i323,2)] ))& + < dptol )) + + !cdp + !1dim + print*,' test_cdp_1dim' + call assert( abs(var(cd1, corrected=.false.) -& + (var(real(cd1), corrected=.false.) +& + var(aimag(cd1), corrected=.false.))) < dptol) + call assert( abs(var(cd1, dim=1, corrected=.false.) -& + (var(real(cd1), dim=1, corrected=.false.) +& + var(aimag(cd1), dim=1, corrected=.false.))) < dptol) + + print*,' test_cdp_1dim_mask' + call assert( isnan(var(cd1, .false., corrected=.false.))) + call assert( isnan(var(cd1, 1, .false., corrected=.false.))) + + print*,' test_cdp_1dim_mask_array' + call assert( abs(var(cd1, aimag(cd1) == 0, corrected=.false.) -& + var(real(cd1), aimag(cd1) == 0, corrected=.false.)) < dptol) + call assert( abs(var(cd1, 1, aimag(cd1) == 0, corrected=.false.) -& + var(real(cd1), 1, aimag(cd1) == 0, corrected=.false.)) < dptol) + call assert( isnan((var(cd1, (real(cd1) == 0 .and. aimag(cd1) == 0),& + corrected=.false.)))) + call assert( abs(var(cd1, (real(cd1) > 1.2 .and. aimag(cd1) == 0),& + corrected=.false.)) < dptol) + + !2dim + cd(:,1) = cd1 + cd(:,2) = cd1*3_sp + cd(:,3) = cd1*1.5_sp + + print*,' test_cdp_2dim' + call assert( abs(var(cd, corrected=.false.) -& + (var(real(cd), corrected=.false.) +& + var(aimag(cd), corrected=.false.))) < dptol) + call assert( all( abs(var(cd, 1, corrected=.false.) -& + (var(real(cd), 1, corrected=.false.) +& + var(aimag(cd), 1, corrected=.false.))) < dptol)) + call assert( all( abs(var(cd, 2, corrected=.false.) -& + (var(real(cd), 2, corrected=.false.) +& + var(aimag(cd), 2, corrected=.false.))) < dptol)) + + print*,' test_cdp_2dim_mask' + call assert( isnan(var(cd, .false., corrected=.false.))) + call assert( any(isnan(var(cd, 1, .false., corrected=.false.)))) + call assert( any(isnan(var(cd, 2, .false., corrected=.false.)))) + + print*,' test_cdp_2dim_mask_array' + call assert( abs(var(cd, aimag(cd) == 0, corrected=.false.) -& + var(real(cd), aimag(cd) == 0, corrected=.false.)) < dptol) + call assert( all( abs( var(cd, 1, aimag(cd) == 0, corrected=.false.) -& + var(real(cd), 1, aimag(cd) == 0, corrected=.false.)) < dptol)) + call assert( any( isnan( var(cd, 2, aimag(cd) == 0, corrected=.false.)))) end program