diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index d4b971ae1..4f3d7e6bb 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -1,5 +1,6 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) +#:set REDRANKS = range(2, MAXRANK + 1) #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_experimental_stats use stdlib_experimental_kinds, only: sp, dp, qp, & @@ -213,9 +214,10 @@ module stdlib_experimental_stats #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_all",rank, t1, k1) - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + ${t1}$, intent(in), optional :: center logical, intent(in), optional :: mask ${t1}$ :: res end function ${RName}$ @@ -225,35 +227,66 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_all",rank, t1, k1, 'dp') - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + real(dp), intent(in), optional :: center logical, intent(in), optional :: mask real(dp) :: res end function ${RName}$ #:endfor #:endfor + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_scalar",rank, t1, k1) + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + ${t1}$, intent(in) :: center + logical, intent(in), optional :: mask + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment",rank, t1, k1) - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in), optional :: mask ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor #:endfor + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_scalar",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + real(dp),intent(in) :: center + logical, intent(in), optional :: mask + real(dp) :: 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("moment",rank, t1, k1, 'dp') - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + real(dp),intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in), optional :: mask real(dp) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ @@ -263,9 +296,10 @@ module stdlib_experimental_stats #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask_all",rank, t1, k1) - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + ${t1}$, intent(in), optional :: center logical, intent(in) :: mask${ranksuffix(rank)}$ ${t1}$ :: res end function ${RName}$ @@ -275,35 +309,66 @@ module stdlib_experimental_stats #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp') - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + real(dp),intent(in), optional :: center logical, intent(in) :: mask${ranksuffix(rank)}$ real(dp) :: res end function ${RName}$ #:endfor #:endfor + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_mask_scalar",rank, t1, k1) + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + ${t1}$, intent(in) :: center + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask",rank, t1, k1) - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in) :: mask${ranksuffix(rank)}$ ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ #:endfor #:endfor + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + real(dp), intent(in) :: center + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: 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("moment_mask",rank, t1, k1, 'dp') - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + real(dp), intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in) :: mask${ranksuffix(rank)}$ real(dp) :: res${reduced_shape('x', rank, 'dim')}$ end function ${RName}$ diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index 95891dc04..c041087a1 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -1,7 +1,7 @@ # Descriptive statistics * [`mean` - mean of array elements](#mean---mean-of-array-elements) -* [`moment` - central moment of array elements](#moment---central-moment-of-array-elements) +* [`moment` - central moments of array elements](#moment---central-moments-of-array-elements) * [`var` - variance of array elements](#var---variance-of-array-elements) ## `mean` - mean of array elements @@ -48,12 +48,15 @@ program demo_mean end program demo_mean ``` -## `moment` - central moment of array elements +## `moment` - central moments of array elements ### Description Returns the _k_-th order central moment 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`. +If a scalar or an array `center` is provided, the function returns the _k_-th order moment about 'center', 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 _k_-th order central moment is defined as : ``` @@ -62,11 +65,17 @@ The _k_-th order central moment is defined as : where n is the number of elements. +The _k_-th order moment about `center` is defined as : + +``` + moment(array) = 1/n sum_i (array(i) - center)^k +``` + ### Syntax -`result = moment(array, order [, mask])` +`result = moment(array, order [, center [, mask]])` -`result = moment(array, order, dim [, mask])` +`result = moment(array, order, dim [, center [, mask]])` ### Arguments @@ -76,6 +85,8 @@ where n is the number of elements. `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`. +`center` (optional): Shall be a scalar of the same type of `result` if `dim` is not provided. If `dim` is provided, `center` shall be a scalar or an array (with a shape similar to that of `array` with dimension `dim` dropped) of the same type of `result`. + `mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`. ### Return value @@ -83,22 +94,24 @@ where n is the number of elements. 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 `real(dp)`. -If `dim` is absent, a scalar with the _k_-th central moment 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 `dim` is absent, a scalar with the _k_-th (central) moment 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 _k_-th central moment of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`. +If `mask` is specified, the result is the _k_-th (central) moment 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_moment - use stdlib_experimental_stats, only: moment + use stdlib_experimental_stats, only: mean, moment implicit none real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ] - print *, moment(x, 2) !returns 2.9167 - print *, moment( reshape(x, [ 2, 3 ] ), 2) !returns 2.9167 - print *, moment( reshape(x, [ 2, 3 ] ), 2, 1) !returns [0.25, 0.25, 0.25] - print *, moment( reshape(x, [ 2, 3 ] ), 2, 1,& - reshape(x, [ 2, 3 ] ) > 3.) !returns [NaN, 0., 0.25] + real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3]) + print *, moment(x, 2) !returns 2.9167 + print *, moment( y, 2) !returns 2.9167 + print *, moment( y, 2, 1) !returns [0.25, 0.25, 0.25] + print *, moment( y, 2, 1, mask = (y > 3.)) !returns [NaN, 0., 0.25] + print *, moment(x, 2, center = 0.) !returns 15.1667 + print *, moment( y, 1, 1, center = 0.) !returns [1.5, 3.5, 5.5] end program demo_moment ``` diff --git a/src/stdlib_experimental_stats_moment.fypp b/src/stdlib_experimental_stats_moment.fypp index 3a6d85dee..1f8afb497 100644 --- a/src/stdlib_experimental_stats_moment.fypp +++ b/src/stdlib_experimental_stats_moment.fypp @@ -1,5 +1,6 @@ #:include "common.fypp" #:set RANKS = range(1, MAXRANK + 1) +#:set REDRANKS = range(2, MAXRANK + 1) #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES submodule (stdlib_experimental_stats) stdlib_experimental_stats_moment @@ -13,14 +14,14 @@ contains #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_all",rank, t1, k1) - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + ${t1}$, intent(in), optional :: center logical, intent(in), optional :: mask ${t1}$ :: res real(${k1}$) :: n - ${t1}$ :: mean if (.not.optval(mask, .true.)) then res = ieee_value(1._${k1}$, ieee_quiet_nan) @@ -28,9 +29,12 @@ contains end if n = size(x, kind = int64) - mean = sum(x) / n - res = sum((x - mean)**order) / n + if (present(center)) then + res = sum((x - center)**order) / n + else + res = sum((x - mean(x))**order) / n + end if end function ${RName}$ #:endfor @@ -40,13 +44,14 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_all",rank, t1, k1, 'dp') - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + real(dp), intent(in), optional :: center logical, intent(in), optional :: mask real(dp) :: res - real(dp) :: n, mean + real(dp) :: n if (.not.optval(mask, .true.)) then res = ieee_value(1._dp, ieee_quiet_nan) @@ -54,9 +59,39 @@ contains end if n = size(x, kind = int64) - mean = sum(real(x, dp)) / n - res = sum((real(x, dp) - mean)**order) / n + if (present(center)) then + res = sum((real(x, dp) - center)**order) / n + else + res = sum((real(x, dp) - mean(x))**order) / n + end if + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_scalar",rank, t1, k1) + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + ${t1}$, intent(in) :: center + logical, intent(in), optional :: mask + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + if (dim >= 1 .and. dim <= ${rank}$) then + res = sum((x - center)**order, dim) / size(x, dim) + else + call error_stop("ERROR (moment): wrong dimension") + end if end function ${RName}$ #:endfor @@ -66,16 +101,17 @@ contains #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment",rank, t1, k1) - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in), optional :: mask ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ integer :: i real(${k1}$) :: n - ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + ${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$ if (.not.optval(mask, .true.)) then res = ieee_value(1._${k1}$, ieee_quiet_nan) @@ -83,15 +119,22 @@ contains end if n = size(x, dim) - mean = sum(x, dim) / n res = 0 select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - do i = 1, size(x, dim) - res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order - end do + if (present(center)) then + do i = 1, size(x, ${fi}$) + res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - center)**order + end do + else + allocate(mean_, source = mean(x, ${fi}$)) + do i = 1, size(x, ${fi}$) + res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order + end do + deallocate(mean_) + end if #:endfor case default call error_stop("ERROR (moment): wrong dimension") @@ -103,19 +146,47 @@ contains #:endfor + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_scalar",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + real(dp),intent(in) :: center + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + if (dim >= 1 .and. dim <= ${rank}$) then + res = sum( (real(x, dp) - center)**order, dim) / size(x, dim) + else + call error_stop("ERROR (moment): wrong dimension") + end if + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment",rank, t1, k1, 'dp') - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + real(dp),intent(in), optional :: center${reduced_shape('x', rank, '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')}$ + real(dp), allocatable :: mean_${ranksuffix(rank-1)}$ if (.not.optval(mask, .true.)) then res = ieee_value(1._dp, ieee_quiet_nan) @@ -123,15 +194,23 @@ contains end if n = size(x, dim) - mean = sum(real(x, dp), dim) / n res = 0 select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - do i = 1, size(x, dim) - res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**order - end do + if (present(center)) then + do i = 1, size(x, ${fi}$) + res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -& + center)**order + end do + else + allocate(mean_, source = mean(x, ${fi}$)) + do i = 1, size(x, ${fi}$) + res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)**order + end do + deallocate(mean_) + end if #:endfor case default call error_stop("ERROR (moment): wrong dimension") @@ -146,19 +225,22 @@ contains #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask_all",rank, t1, k1) - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + ${t1}$, intent(in), optional :: center logical, intent(in) :: mask${ranksuffix(rank)}$ ${t1}$ :: res real(${k1}$) :: n - ${t1}$ :: mean n = count(mask, kind = int64) - mean = sum(x, mask) / n - res = sum((x - mean)**order, mask) / n + if (present(center)) then + res = sum((x - center)**order, mask) / n + else + res = sum((x - mean(x, mask))**order, mask) / n + end if end function ${RName}$ #:endfor @@ -168,18 +250,44 @@ contains #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp') - module function ${RName}$(x, order, mask) result(res) + module function ${RName}$(x, order, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order + real(dp),intent(in), optional :: center logical, intent(in) :: mask${ranksuffix(rank)}$ real(dp) :: res - real(dp) :: n, mean + real(dp) :: n n = count(mask, kind = int64) - mean = sum(real(x, dp), mask) / n - res = sum((real(x, dp) - mean)**order, mask) / n + if (present(center)) then + res = sum((real(x, dp) - center)**order, mask) / n + else + res = sum((real(x, dp) - mean(x,mask))**order, mask) / n + end if + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_mask_scalar",rank, t1, k1) + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + ${t1}$, intent(in) :: center + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + + if (dim >= 1 .and. dim <= ${rank}$) then + res = sum((x - center)**order, dim, mask) / count(mask, dim) + else + call error_stop("ERROR (moment): wrong dimension") + end if end function ${RName}$ #:endfor @@ -189,33 +297,48 @@ contains #:for k1, t1 in RC_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask",rank, t1, k1) - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$ logical, intent(in) :: mask${ranksuffix(rank)}$ ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ integer :: i real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$ - ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + ${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$ n = count(mask, dim) - mean = sum(x, dim, mask) / n res = 0 select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - do i = 1, size(x, dim) - res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order,& - #:if t1[0] == 'r' - 0._${k1}$,& - #:else - cmplx(0,0,kind=${k1}$),& - #:endif - mask${select_subarray(rank, [(fi, 'i')])}$) - end do + if (present(center)) then + do i = 1, size(x, ${fi}$) + res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ -& + center)**order,& + #:if t1[0] == 'r' + 0._${k1}$,& + #:else + cmplx(0,0,kind=${k1}$),& + #:endif + mask${select_subarray(rank, [(fi, 'i')])}$) + end do + else + allocate(mean_, source = mean(x, ${fi}$, mask)) + do i = 1, size(x, ${fi}$) + res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order,& + #:if t1[0] == 'r' + 0._${k1}$,& + #:else + cmplx(0,0,kind=${k1}$),& + #:endif + mask${select_subarray(rank, [(fi, 'i')])}$) + end do + deallocate(mean_) + end if #:endfor case default call error_stop("ERROR (moment): wrong dimension") @@ -227,31 +350,64 @@ contains #:endfor + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in REDRANKS + #:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, center, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + real(dp), intent(in) :: center + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + if (dim >= 1 .and. dim <= ${rank}$) then + res = sum(( real(x, dp) - center)**order, dim, mask) / count(mask, dim) + else + call error_stop("ERROR (moment): wrong dimension") + end if + + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES #:for rank in RANKS #:set RName = rname("moment_mask",rank, t1, k1, 'dp') - module function ${RName}$(x, order, dim, mask) result(res) + module function ${RName}$(x, order, dim, center, mask) result(res) ${t1}$, intent(in) :: x${ranksuffix(rank)}$ integer, intent(in) :: order integer, intent(in) :: dim + real(dp), intent(in), optional :: center${reduced_shape('x', rank, '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')}$ + real(dp), allocatable :: mean_${ranksuffix(rank-1)}$ n = count(mask, dim) - mean = sum(real(x, dp), dim, mask) / n res = 0 select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - do i = 1, size(x, dim) - res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**order,& - 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) - end do + if (present(center)) then + do i = 1, size(x, ${fi}$) + res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -& + center)**order,& + 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) + end do + else + allocate(mean_, source = mean(x, ${fi}$, mask)) + do i = 1, size(x, ${fi}$) + res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)& + **order,& + 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) + end do + deallocate(mean_) + end if #:endfor case default call error_stop("ERROR (moment): wrong dimension") diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index aa0a9e1ba..b79e39bf6 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,5 +1,6 @@ ADDTEST(mean) ADDTEST(moment) +ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) diff --git a/src/tests/stats/test_moment.f90 b/src/tests/stats/test_moment.f90 index 8eea4529d..fd00d51a1 100644 --- a/src/tests/stats/test_moment.f90 +++ b/src/tests/stats/test_moment.f90 @@ -49,12 +49,12 @@ subroutine test_sp(x1, x2) call check( abs(moment(x1, order, dim=1)) < sptol) print*,' test_sp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_sp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5)) < sptol) - call check( abs(moment(x1, order, 1, x1 < 5)) < sptol) + call check( abs(moment(x1, order, mask = (x1 < 5))) < sptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5))) < sptol) !2dim print*,' test_sp_2dim', order @@ -63,14 +63,14 @@ subroutine test_sp(x1, x2) call check( all( abs( moment(x2, order, 2)) < sptol)) print*,' test_sp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_sp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)) < sptol) - call check( all( abs( moment(x2, order, 1, x2 < 11)) < sptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11)) < sptol)) + call check( abs(moment(x2, order, mask = (x2 < 11))) < sptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11))) < sptol)) + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11))) < sptol)) !3dim allocate(x3(size(x2,1),size(x2,2),3)) @@ -85,16 +85,16 @@ subroutine test_sp(x1, x2) call check( all( abs( moment(x3, order, 3)) < sptol)) print*,' test_sp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_sp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) ) < sptol) - call check( all( abs( moment(x3, order, 1, x3 < 45)) < sptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45)) < sptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45)) < sptol )) + call check( abs(moment(x3, order, mask = (x3 < 11)) ) < sptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45))) < sptol )) + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45))) < sptol )) + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45))) < sptol )) order = 2 @@ -105,29 +105,30 @@ subroutine test_sp(x1, x2) call check( abs(moment(x1, order, dim=1) - 2._sp) < sptol) print*,' test_sp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_sp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5) - 1.25_sp) < sptol) - call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_sp) < sptol) + call check( abs(moment(x1, order, mask = (x1 < 5)) - 1.25_sp) < sptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5)) - 1.25_sp) < sptol) !2dim print*,' test_sp_2dim', order call check( abs(moment(x2, order) - 107.25_sp/9.) < sptol) call check( all( abs( moment(x2, order, 1) - [5._sp, 5._sp, 1.25_sp]) < sptol)) - call check( all( abs( moment(x2, order, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol)) + call check( all( abs( moment(x2, order, 2) -& + [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol)) print*,' test_sp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_sp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)- 2.75_sp*3.) < sptol) - call check( all( abs( moment(x2, order, 1, x2 < 11) -& + call check( abs(moment(x2, order, mask = (x2 < 11))- 2.75_sp*3.) < sptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11)) -& [5._sp, 5._sp, 0.25_sp]) < sptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11) -& + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11)) -& [19._sp*2./3., 43._sp/9.*2., 0.25_sp , 0.25_sp]) < sptol)) !3dim @@ -154,24 +155,25 @@ subroutine test_sp(x1, x2) < sptol)) print*,' test_sp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_sp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < sptol) - call check( all( abs( moment(x3, order, 1, x3 < 45) -& + call check( abs(moment(x3, order, mask = (x3 < 11)) -& + 7.7370242214532876_dp ) < sptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45)) -& reshape([5._sp, 5._sp, 1.25_sp, 20._sp, 20._sp, 5._sp,& 80._sp, 80._sp, 32._sp/3.],& [size(x3, 2), size(x3, 3)])) < sptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45) -& + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45)) -& reshape([ 38._sp/3., 86._sp/9., 62._sp/9., 14._sp/3., 152._sp/3.,& 344._sp/9., 248._sp/9., 168._sp/9., 1824._sp/9.,& 1376._sp/9., 992._sp/9., 4._sp& ],& [size(x3, 1), size(x3, 3)])) < sptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45) -& + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45)) -& reshape([14._sp/9., 14._sp, 350._sp/9., 686._sp/9., 56._sp/9.,& 224._sp/9., 56._sp, 896._sp/9., 126._sp, 1400._sp/9.,& 1694._sp/9., 36._sp& @@ -194,12 +196,12 @@ subroutine test_dp(x1, x2) call check( abs(moment(x1, order, dim=1)) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5)) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5)) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5))) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5))) < dptol) !2dim print*,' test_dp_2dim', order @@ -208,14 +210,14 @@ subroutine test_dp(x1, x2) call check( all( abs( moment(x2, order, 2)) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + call check( abs(moment(x2, order, mask = (x2 < 11))) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11))) < dptol)) + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11))) < dptol)) !3dim allocate(x3(size(x2,1),size(x2,2),3)) @@ -230,16 +232,16 @@ subroutine test_dp(x1, x2) call check( all( abs( moment(x3, order, 3)) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + call check( abs(moment(x3, order, mask = (x3 < 11)) ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45))) < dptol )) order = 2 @@ -250,12 +252,12 @@ subroutine test_dp(x1, x2) call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5)) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5)) - 1.25_dp) < dptol) !2dim print*,' test_dp_2dim', order @@ -265,15 +267,15 @@ subroutine test_dp(x1, x2) [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11) -& + call check( abs(moment(x2, order, mask = (x2 < 11))- 2.75_dp*3.) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11)) -& [5._dp, 5._dp, 0.25_dp]) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11) -& + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11)) -& [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) !3dim @@ -300,24 +302,25 @@ subroutine test_dp(x1, x2) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45) -& + call check( abs(moment(x3, order, mask = (x3 < 11)) -& + 7.7370242214532876_dp ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45)) -& reshape([5._dp, 5._dp, 1.25_dp, 20._dp, 20._dp, 5._dp,& 80._dp, 80._dp, 32._dp/3.],& [size(x3, 2), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45) -& + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45)) -& 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(x3, 1), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45) -& + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45)) -& 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& @@ -340,12 +343,12 @@ subroutine test_int32(x1, x2) call check( abs(moment(x1, order, dim=1)) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5)) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5)) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5))) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5))) < dptol) !2dim print*,' test_dp_2dim', order @@ -354,14 +357,14 @@ subroutine test_int32(x1, x2) call check( all( abs( moment(x2, order, 2)) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + call check( abs(moment(x2, order, mask = (x2 < 11))) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11))) < dptol)) + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11))) < dptol)) !3dim allocate(x3(size(x2,1),size(x2,2),3)) @@ -376,16 +379,16 @@ subroutine test_int32(x1, x2) call check( all( abs( moment(x3, order, 3)) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + call check( abs(moment(x3, order, mask = (x3 < 11)) ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45))) < dptol )) order = 2 @@ -396,12 +399,12 @@ subroutine test_int32(x1, x2) call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5)) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5)) - 1.25_dp) < dptol) !2dim print*,' test_dp_2dim', order @@ -411,15 +414,15 @@ subroutine test_int32(x1, x2) [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11) -& + call check( abs(moment(x2, order, mask = (x2 < 11))- 2.75_dp*3.) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11)) -& [5._dp, 5._dp, 0.25_dp]) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11) -& + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11)) -& [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) !3dim @@ -446,24 +449,25 @@ subroutine test_int32(x1, x2) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45) -& + call check( abs(moment(x3, order, mask = (x3 < 11)) -& + 7.7370242214532876_dp ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45)) -& reshape([5._dp, 5._dp, 1.25_dp, 20._dp, 20._dp, 5._dp,& 80._dp, 80._dp, 32._dp/3.],& [size(x3, 2), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45) -& + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45)) -& 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(x3, 1), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45) -& + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45)) -& 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& @@ -486,12 +490,12 @@ subroutine test_int64(x1, x2) call check( abs(moment(x1, order, dim=1)) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5)) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5)) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5))) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5))) < dptol) !2dim print*,' test_dp_2dim', order @@ -500,14 +504,14 @@ subroutine test_int64(x1, x2) call check( all( abs( moment(x2, order, 2)) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + call check( abs(moment(x2, order, mask = (x2 < 11))) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11))) < dptol)) + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11))) < dptol)) !3dim allocate(x3(size(x2,1),size(x2,2),3)) @@ -522,16 +526,16 @@ subroutine test_int64(x1, x2) call check( all( abs( moment(x3, order, 3)) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + call check( abs(moment(x3, order, mask = (x3 < 11)) ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45))) < dptol )) + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45))) < dptol )) order = 2 @@ -542,12 +546,12 @@ subroutine test_int64(x1, x2) call check( abs(moment(x1, order, dim=1) - 2._dp) < dptol) print*,' test_dp_1dim_mask', order - call check( ieee_is_nan(moment(x1, order, .false.))) - call check( ieee_is_nan(moment(x1, order, 1, .false.))) + call check( ieee_is_nan(moment(x1, order, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, 1, mask = .false.))) print*,' test_dp_1dim_mask_array', order - call check( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) - call check( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, mask = (x1 < 5)) - 1.25_dp) < dptol) + call check( abs(moment(x1, order, 1, mask = (x1 < 5)) - 1.25_dp) < dptol) !2dim print*,' test_dp_2dim', order @@ -557,15 +561,15 @@ subroutine test_int64(x1, x2) [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) print*,' test_dp_2dim_mask', order - call check( ieee_is_nan(moment(x2, order, .false.))) - call check( any(ieee_is_nan(moment(x2, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x2, order, 2, .false.)))) + call check( ieee_is_nan(moment(x2, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, 2, mask = .false.)))) print*,' test_dp_2dim_mask_array', order - call check( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) - call check( all( abs( moment(x2, order, 1, x2 < 11) -& + call check( abs(moment(x2, order, mask = (x2 < 11))- 2.75_dp*3.) < dptol) + call check( all( abs( moment(x2, order, 1, mask = (x2 < 11)) -& [5._dp, 5._dp, 0.25_dp]) < dptol)) - call check( all( abs( moment(x2, order, 2, x2 < 11) -& + call check( all( abs( moment(x2, order, 2, mask = (x2 < 11)) -& [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) !3dim @@ -592,24 +596,25 @@ subroutine test_int64(x1, x2) < dptol)) print*,' test_dp_3dim_mask', order - call check( ieee_is_nan(moment(x3, order, .false.))) - call check( any(ieee_is_nan(moment(x3, order, 1, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 2, .false.)))) - call check( any(ieee_is_nan(moment(x3, order, 3, .false.)))) + call check( ieee_is_nan(moment(x3, order, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, 1, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 2, mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, 3, mask = .false.)))) print*,' test_dp_3dim_mask_array', order - call check( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) - call check( all( abs( moment(x3, order, 1, x3 < 45) -& + call check( abs(moment(x3, order, mask = (x3 < 11)) -& + 7.7370242214532876_dp ) < dptol) + call check( all( abs( moment(x3, order, 1, mask = (x3 < 45)) -& reshape([5._dp, 5._dp, 1.25_dp, 20._dp, 20._dp, 5._dp,& 80._dp, 80._dp, 32._dp/3.],& [size(x3, 2), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 2, x3 < 45) -& + call check( all( abs( moment(x3, order, 2, mask = (x3 < 45)) -& 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(x3, 1), size(x3, 3)])) < dptol )) - call check( all( abs( moment(x3, order, 3, x3 < 45) -& + call check( all( abs( moment(x3, order, 3, mask = (x3 < 45)) -& 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& @@ -622,7 +627,6 @@ subroutine test_csp(x1, x2) complex(sp), intent(in) :: x1(:), x2(:, :) integer :: order - complex(sp), allocatable :: x3(:, :, :) order = 1 @@ -632,12 +636,12 @@ subroutine test_csp(x1, x2) call check( abs(moment(x1, order, dim=1)) < sptol) print*,' test_sp_1dim_mask', order - call check( ieee_is_nan(abs(moment(x1, order, .false.)))) - call check( ieee_is_nan(abs(moment(x1, order, 1, .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, mask = .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, 1, mask = .false.)))) print*,' test_sp_1dim_mask_array', order - call check( abs(moment(x1, order, aimag(x1) == 0)) < sptol) - call check( abs(moment(x1, order, 1, aimag(x1) == 0)) < sptol) + call check( abs(moment(x1, order, mask = (aimag(x1) == 0))) < sptol) + call check( abs(moment(x1, order, 1, mask = (aimag(x1) == 0))) < sptol) !2dim print*,' test_sp_2dim', order @@ -646,14 +650,15 @@ subroutine test_csp(x1, x2) call check( all( abs( moment(x2, order, 2)) < sptol)) print*,' test_sp_2dim_mask', order - call check( ieee_is_nan(abs(moment(x2, order, .false.)))) - call check( any(ieee_is_nan(abs(moment(x2, order, 1, .false.))))) - call check( any(ieee_is_nan(abs(moment(x2, order, 2, .false.))))) + call check( ieee_is_nan(abs(moment(x2, order, mask = .false.)))) + call check( any(ieee_is_nan(abs(moment(x2, order, 1, mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, 2, mask = .false.))))) print*,' test_sp_2dim_mask_array', order - call check( abs(moment(x2, order, aimag(x2) == 0)) < sptol) - call check( all( abs( moment(x2, order, 1, aimag(x2) == 0)) < sptol)) - call check( any(ieee_is_nan( abs( moment(x2, order, 2, aimag(x2) == 0))))) + call check( abs(moment(x2, order, mask = (aimag(x2) == 0))) < sptol) + call check( all( abs( moment(x2, order, 1, mask = (aimag(x2) == 0))) < sptol)) + call check( any(ieee_is_nan( abs( moment(x2, order, 2,& + mask = (aimag(x2) == 0)))))) order = 2 @@ -664,13 +669,13 @@ subroutine test_csp(x1, x2) (-6.459422410E-02,-0.556084037)) < sptol) print*,' test_sp_1dim_mask', order - call check( ieee_is_nan(abs(moment(x1, order, .false.)))) - call check( ieee_is_nan(abs(moment(x1, order, 1, .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, mask = .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, 1, mask = .false.)))) print*,' test_sp_1dim_mask_array', order - call check( abs(moment(x1, order, aimag(x1) == 0) -& + call check( abs(moment(x1, order, mask = (aimag(x1) == 0)) -& (8.969944715E-02,0.00000000)) < sptol) - call check( abs(moment(x1, order, 1, aimag(x1) == 0) -& + call check( abs(moment(x1, order, 1, mask = (aimag(x1) == 0)) -& (8.969944715E-02,0.00000000)) < sptol) !2dim @@ -690,14 +695,14 @@ subroutine test_csp(x1, x2) ) < sptol)) print*,' test_sp_2dim_mask', order - call check( ieee_is_nan(abs(moment(x2, order, .false.)))) - call check( any(ieee_is_nan(abs(moment(x2, order, 1, .false.))))) - call check( any(ieee_is_nan(abs(moment(x2, order, 2, .false.))))) + call check( ieee_is_nan(abs(moment(x2, order, mask = .false.)))) + call check( any(ieee_is_nan(abs(moment(x2, order, 1, mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, 2, mask = .false.))))) print*,' test_sp_2dim_mask_array', order - call check( abs(moment(x2, order, aimag(x2) == 0)-& + call check( abs(moment(x2, order, mask = (aimag(x2) == 0))-& (1.08109438,0.00000000)) < sptol) - call check( all( abs( moment(x2, order, 1, aimag(x2)==0) -& + call check( all( abs( moment(x2, order, 1, mask = (aimag(x2)==0)) -& [(8.969944715E-02,0.00000000),& (0.807295084,0.00000000),& (0.201823771,0.00000000)]& diff --git a/src/tests/stats/test_rawmoment.f90 b/src/tests/stats/test_rawmoment.f90 new file mode 100644 index 000000000..58cdd8d3b --- /dev/null +++ b/src/tests/stats/test_rawmoment.f90 @@ -0,0 +1,624 @@ +program test_rawmoment + use stdlib_experimental_error, only: check + use stdlib_experimental_kinds, only: sp, dp, int32, int64 + use stdlib_experimental_stats, only: mean, moment + use,intrinsic :: ieee_arithmetic, only : ieee_is_nan + implicit none + + + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 1000 * epsilon(1._dp) + + real(dp) :: d1(5) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_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(sp) :: cs(5,3) + + + call test_sp(real(d1,sp), real(d,sp)) + call test_int32(int(d1, int32), int(d, int32)) + + cs(:,1) = cs1 + cs(:,2) = cs1*3_sp + cs(:,3) = cs1*1.5_sp + call test_csp(cs1, cs) + + +contains + subroutine test_sp(x1, x2) + real(sp), intent(in) :: x1(:), x2(:, :) + + integer :: order + real(sp), allocatable :: x3(:, :, :) + real(sp), allocatable :: zero2_1(:), zero2_2(:) + real(sp), allocatable :: zero3_1(:,:), zero3_2(:,:), zero3_3(:,:) + + allocate(zero2_1(size(x2, 2)), zero2_2(size(x2, 1))) + zero2_1 = 0 + zero2_2 = 0 + + order = 1 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = 0.) - mean(x1)) < sptol) + call check( abs(moment(x1, order, 1, center = 0.) - mean(x1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(moment(x1, order, center = 0., mask = .false.))) + call check( ieee_is_nan(moment(x1, order, dim = 1, center = 0., mask = .false.))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = 0., mask = (x1 < 5)) -& + mean(x1, mask = (x1 < 5)) ) < sptol) + call check( abs(moment(x1, order, dim = 1, center = 0., mask = (x1 < 5)) -& + mean(x1, dim = 1, mask = (x1 < 5))) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = 0.) - mean(x2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0.) -& + mean(x2, dim = 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0.) -& + mean(x2, dim = 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2, dim = 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) -& + mean(x2, dim = 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(moment(x2, order, center = 0., mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, dim = 1, center = zero2_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 2, center = zero2_2,& + mask = .false.)))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = 0., mask = (x2 < 11)) -& + mean(x2, x2 < 11)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0.,& + mask = (x2 < 11)) -& + mean(x2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0.,& + mask = (x2 < 11)) -& + mean(x2, 2, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (x2 < 11)) -& + mean(x2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2,& + mask = (x2 < 11)) -& + mean(x2, 2, x2 < 11)) < sptol)) + + !3dim + allocate(x3(size(x2,1),size(x2,2),3)) + x3(:,:,1)=x2; + x3(:,:,2)=x2*2; + x3(:,:,3)=x2*4; + + allocate(zero3_1(size(x3, 2), size(x3, 3))& + ,zero3_2(size(x3, 1), size(x3, 3))& + ,zero3_3(size(x3, 1), size(x3, 2))) + zero3_1 = 0 + zero3_2 = 0 + zero3_3 = 0 + + print*,' test_sp_3dim', order + call check( abs(moment(x3, order, center = 0.) - mean(x3)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0._sp) -& + mean(x3, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = 0._sp) -& + mean(x3, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = 0._sp) -& + mean(x3, 3)) < sptol)) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1) -& + mean(x3, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2) -& + mean(x3, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3) -& + mean(x3, 3)) < sptol)) + + print*,' test_sp_3dim_mask', order + call check( ieee_is_nan(moment(x3, order, center = 0., mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, dim = 1, center = zero3_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 2, center = zero3_2,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 3, center = zero3_3,& + mask = .false.)))) + + print*,' test_sp_3dim_mask_array', order + call check( abs(moment(x3, order, center = 0., mask = (x3 < 11)) -& + mean(x3, x3 < 11)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0.,& + mask = (x3 < 45)) -& + mean(x3, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = 0.,& + mask = (x3 < 45)) -& + mean(x3, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = 0.,& + mask = (x3 < 45)) -& + mean(x3, 3, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1,& + mask = (x3 < 45)) -& + mean(x3, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2,& + mask = (x3 < 45)) -& + mean(x3, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3,& + mask = (x3 < 45)) -& + mean(x3, 3, x3 < 45)) < sptol )) + + + order = 2 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = 0.) - mean(x1**2)) < sptol) + call check( abs(moment(x1, order, 1, center = 0.) - mean(x1**2, 1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(moment(x1, order, center = 0., mask = .false.))) + call check( ieee_is_nan(moment(x1, order, dim = 1, center = 0., mask = .false.))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = 0., mask = (x1 < 5)) -& + mean(x1**2, x1 < 5)) < sptol) + call check( abs(moment(x1, order, dim = 1, center = 0., mask = (x1 < 5)) -& + mean(x1**2, 1, x1 < 5)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = 0.) - mean(x2**2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0.) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0.) - & + mean(x2**2, 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) - & + mean(x2**2, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(moment(x2, order, center = 0., mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, dim = 1, center = zero2_1, & + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 2, center = zero2_2, & + mask = .false.)))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = 0., mask = (x2 < 11)) -& + mean(x2**2, x2 < 11)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0.,& + mask = (x2 < 11)) -& + mean(x2**2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0.,& + mask = (x2 < 11)) -& + mean(x2**2, 2, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (x2 < 11)) -& + mean(x2**2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2,& + mask = (x2 < 11)) -& + mean(x2**2, 2, x2 < 11)) < sptol)) + + !3dim + print*,' test_sp_3dim', order + call check( abs(moment(x3, order, center = 0.) - mean(x3**2)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0.) -& + mean(x3**2, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = 0.) -& + mean(x3**2, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = 0.) -& + mean(x3**2, 3)) < sptol)) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1) -& + mean(x3**2, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2) -& + mean(x3**2, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3) -& + mean(x3**2, 3)) < sptol)) + + print*,' test_sp_3dim_mask', order + call check( ieee_is_nan(moment(x3, order, center = 0., mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, dim = 1, center = zero3_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 2, center = zero3_2,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 3, center = zero3_3,& + mask = .false.)))) + + print*,' test_sp_3dim_mask_array', order + call check( abs(moment(x3, order, center = 0., mask = (x3 < 11)) -& + mean(x3**2, x3 < 11)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0.,& + mask = (x3 < 45)) -& + mean(x3**2, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = 0.,& + mask = (x3 < 45)) -& + mean(x3**2, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = 0.,& + mask = (x3 < 45)) -& + mean(x3**2, 3, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1,& + mask = (x3 < 45)) -& + mean(x3**2, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2,& + mask = (x3 < 45)) -& + mean(x3**2, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3,& + mask = (x3 < 45)) -& + mean(x3**2, 3, x3 < 45)) < sptol )) + + end subroutine + + subroutine test_int32(x1, x2) + integer(int32), intent(in) :: x1(:), x2(:, :) + + integer :: order + integer(int32), allocatable :: x3(:, :, :) + real(dp), allocatable :: zero2_1(:), zero2_2(:) + real(dp), allocatable :: zero3_1(:,:), zero3_2(:,:), zero3_3(:,:) + + allocate(zero2_1(size(x2, 2)), zero2_2(size(x2, 1))) + zero2_1 = 0 + zero2_2 = 0 + + order = 1 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = 0._dp) - mean(x1)) < sptol) + call check( abs(moment(x1, order, 1, center = 0._dp) - mean(x1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(moment(x1, order, center = 0._dp, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, dim = 1, center = 0._dp,& + mask = .false.))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = 0._dp, mask = (x1 < 5)) -& + mean(x1, mask = (x1 < 5)) ) < sptol) + call check( abs(moment(x1, order, dim = 1, center = 0._dp, mask = (x1 < 5)) -& + mean(x1, dim = 1, mask = (x1 < 5))) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = 0._dp) - mean(x2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0._dp) -& + mean(x2, dim = 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0._dp) -& + mean(x2, dim = 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2, dim = 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) -& + mean(x2, dim = 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(moment(x2, order, center = 0._dp, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, dim = 1, center = zero2_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 2, center = zero2_2,& + mask = .false.)))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = 0._dp, mask = (x2 < 11)) -& + mean(x2, x2 < 11)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0._dp,& + mask = (x2 < 11)) -& + mean(x2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0._dp,& + mask = (x2 < 11)) -& + mean(x2, 2, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (x2 < 11)) -& + mean(x2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2,& + mask = (x2 < 11)) -& + mean(x2, 2, x2 < 11)) < sptol)) + + !3dim + allocate(x3(size(x2,1),size(x2,2),3)) + x3(:,:,1)=x2; + x3(:,:,2)=x2*2; + x3(:,:,3)=x2*4; + + allocate(zero3_1(size(x3, 2), size(x3, 3))& + ,zero3_2(size(x3, 1), size(x3, 3))& + ,zero3_3(size(x3, 1), size(x3, 2))) + zero3_1 = 0 + zero3_2 = 0 + zero3_3 = 0 + + print*,' test_sp_3dim', order + call check( abs(moment(x3, order, center = 0._dp) - mean(x3)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0._dp) -& + mean(x3, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = 0._dp) -& + mean(x3, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = 0._dp) -& + mean(x3, 3)) < sptol)) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1) -& + mean(x3, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2) -& + mean(x3, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3) -& + mean(x3, 3)) < sptol)) + + print*,' test_sp_3dim_mask', order + call check( ieee_is_nan(moment(x3, order, center = 0._dp, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, dim = 1, center = zero3_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 2, center = zero3_2,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 3, center = zero3_3,& + mask = .false.)))) + + print*,' test_sp_3dim_mask_array', order + call check( abs(moment(x3, order, center = 0._dp, mask = (x3 < 11)) -& + mean(x3, x3 < 11)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3, 3, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1,& + mask = (x3 < 45)) -& + mean(x3, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2,& + mask = (x3 < 45)) -& + mean(x3, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3,& + mask = (x3 < 45)) -& + mean(x3, 3, x3 < 45)) < sptol )) + + + order = 2 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = 0._dp) - mean(x1**2)) < sptol) + call check( abs(moment(x1, order, 1, center = 0._dp) - mean(x1**2, 1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(moment(x1, order, center = 0._dp, mask = .false.))) + call check( ieee_is_nan(moment(x1, order, dim = 1, center = 0._dp,& + mask = .false.))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = 0._dp, mask = (x1 < 5)) -& + mean(x1**2, x1 < 5)) < sptol) + call check( abs(moment(x1, order, dim = 1, center = 0._dp, mask = (x1 < 5)) -& + mean(x1**2, 1, x1 < 5)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = 0._dp) - mean(x2**2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0._dp) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0._dp) -& + mean(x2**2, 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) -& + mean(x2**2, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(moment(x2, order, center = 0._dp, mask = .false.))) + call check( any(ieee_is_nan(moment(x2, order, dim = 1, center = 0._dp,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 2, center = 0._dp,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 1, center = zero2_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x2, order, dim = 2, center = zero2_2,& + mask = .false.)))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = 0._dp, mask = (x2 < 11)) -& + mean(x2**2, x2 < 11)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = 0._dp,& + mask = (x2 < 11)) -& + mean(x2**2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = 0._dp,& + mask = (x2 < 11)) -& + mean(x2**2, 2, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (x2 < 11)) -& + mean(x2**2, 1, x2 < 11)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2,& + mask = (x2 < 11)) -& + mean(x2**2, 2, x2 < 11)) < sptol)) + + !3dim + print*,' test_sp_3dim', order + call check( abs(moment(x3, order, center = 0._dp) - mean(x3**2)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0._dp) -& + mean(x3**2, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = 0._dp) -& + mean(x3**2, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = 0._dp) -& + mean(x3**2, 3)) < sptol)) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1) -& + mean(x3**2, 1)) < sptol)) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2) -& + mean(x3**2, 2)) < sptol)) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3) -& + mean(x3**2, 3)) < sptol)) + + print*,' test_sp_3dim_mask', order + call check( ieee_is_nan(moment(x3, order, center = 0._dp, mask = .false.))) + call check( any(ieee_is_nan(moment(x3, order, dim = 1, center = 0._dp,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 2, center = 0._dp,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 3, center = 0._dp,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 1, center = zero3_1,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 2, center = zero3_2,& + mask = .false.)))) + call check( any(ieee_is_nan(moment(x3, order, dim = 3, center = zero3_3,& + mask = .false.)))) + + print*,' test_sp_3dim_mask_array', order + call check( abs(moment(x3, order, center = 0._dp, mask = (x3 < 11)) -& + mean(x3**2, x3 < 11)) < sptol) + call check( all( abs( moment(x3, order, dim = 1, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3**2, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3**2, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = 0._dp,& + mask = (x3 < 45)) -& + mean(x3**2, 3, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 1, center = zero3_1,& + mask = (x3 < 45)) -& + mean(x3**2, 1, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 2, center = zero3_2,& + mask = (x3 < 45)) -& + mean(x3**2, 2, x3 < 45)) < sptol )) + call check( all( abs( moment(x3, order, dim = 3, center = zero3_3,& + mask = (x3 < 45)) -& + mean(x3**2, 3, x3 < 45)) < sptol )) + + end subroutine + + subroutine test_csp(x1, x2) + complex(sp), intent(in) :: x1(:), x2(:, :) + + integer :: order + complex(sp), allocatable :: zero2_1(:), zero2_2(:) + + allocate(zero2_1(size(x2, 2)), zero2_2(size(x2, 1))) + zero2_1 = (0., 0.) + zero2_2 = (0., 0.) + + order = 1 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = (0., 0.)) - mean(x1)) < sptol) + call check( abs(moment(x1, order, 1, center = (0., 0.)) - mean(x1, 1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(abs(moment(x1, order, center = (0., 0.),& + mask = .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, dim = 1, center = (0., 0.),& + mask = .false.)))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = (0., 0.), mask = (aimag(x1) == 0)) -& + mean(x1, aimag(x1) == 0)) < sptol) + call check( abs(moment(x1, order, dim = 1, center = (0., 0.),& + mask = (aimag(x1) == 0)) -& + mean(x1, 1, aimag(x1) == 0)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = (0., 0.)) - mean(x2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = (0., 0.)) -& + mean(x2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = (0., 0.)) -& + mean(x2, 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) -& + mean(x2, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(abs(moment(x2, order, center = (0., 0.),& + mask = .false.)))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 1, center = (0., 0.),& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 2, center = (0., 0.),& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 1, center = zero2_1,& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 2, center = zero2_2,& + mask = .false.))))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = (0., 0.), mask = (aimag(x2) == 0)) -& + mean(x2, aimag(x2) == 0)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = (0., 0.),& + mask = (aimag(x2) == 0)) -& + mean(x2, 1, aimag(x2) == 0)) < sptol)) + call check( any(ieee_is_nan( abs( moment(x2, order,& + dim = 2, center = (0., 0.), mask = (aimag(x2) == 0)) -& + mean(x2, 2, aimag(x2) == 0))))) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (aimag(x2) == 0)) -& + mean(x2, 1, aimag(x2) == 0)) < sptol)) + call check( any(ieee_is_nan( abs( moment(x2, order,& + dim = 2, center = zero2_2, mask = (aimag(x2) == 0)) -& + mean(x2, 2, aimag(x2) == 0))))) + + order = 2 + + !1dim + print*,' test_sp_1dim', order + call check( abs(moment(x1, order, center = (0., 0.)) - mean(x1**2)) < sptol) + call check( abs(moment(x1, order, 1, center = (0., 0.)) -& + mean(x1**2, 1)) < sptol) + + print*,' test_sp_1dim_mask', order + call check( ieee_is_nan(abs(moment(x1, order, center = (0., 0.),& + mask = .false.)))) + call check( ieee_is_nan(abs(moment(x1, order, dim = 1, center = (0., 0.),& + mask = .false.)))) + + print*,' test_sp_1dim_mask_array', order + call check( abs(moment(x1, order, center = (0., 0.), mask = (aimag(x1) == 0)) -& + mean(x1**2, aimag(x1) == 0)) < sptol) + call check( abs(moment(x1, order, dim = 1, center = (0., 0.),& + mask = (aimag(x1) == 0)) -& + mean(x1**2, 1, aimag(x1) == 0)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call check( abs(moment(x2, order, center = (0., 0.)) - mean(x2**2)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = (0., 0.)) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = (0., 0.)) -& + mean(x2**2, 2)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1) -& + mean(x2**2, 1)) < sptol)) + call check( all( abs( moment(x2, order, dim = 2, center = zero2_2) -& + mean(x2**2, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call check( ieee_is_nan(abs(moment(x2, order, center = (0., 0.),& + mask = .false.)))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 1, center = (0., 0.),& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 2, center = (0., 0.),& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 1, center = zero2_1,& + mask = .false.))))) + call check( any(ieee_is_nan(abs(moment(x2, order, dim = 2, center = zero2_2,& + mask = .false.))))) + + print*,' test_sp_2dim_mask_array', order + call check( abs(moment(x2, order, center = (0., 0.), mask = (aimag(x2) == 0)) -& + mean(x2**2, aimag(x2) == 0)) < sptol) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (aimag(x2)==0)) -& + mean(x2**2, 1, aimag(x2)==0)) < sptol)) + call check( all( abs( moment(x2, order, dim = 1, center = zero2_1,& + mask = (aimag(x2)==0)) -& + mean(x2**2, 1, aimag(x2)==0)) < sptol)) + + end subroutine +end program