From de371129f3a9af9b102e0f313d321574a0505c32 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Sun, 1 Mar 2020 20:22:53 +0100 Subject: [PATCH 1/3] Squashed commit of the following: commit c9f73b52d98600b1783c2f4a662ff8fba88cb525 Author: Vandenplas, Jeremie Date: Sun Mar 1 20:20:58 2020 +0100 moment_dev commit df4fedbf231830297e64c719025df70ca96e601c Author: Vandenplas, Jeremie Date: Sat Feb 29 22:26:53 2020 +0100 moment_dev: addition of TOCs commit 20a76b82de37a340c66b48e832afec68277ac05e Author: Vandenplas, Jeremie Date: Sat Feb 29 22:19:03 2020 +0100 moment_dev: addition of spec commit 7c33e00cbe1d1efb91f6fb59ad27043582ee6f00 Author: Vandenplas, Jeremie Date: Sat Feb 29 21:07:24 2020 +0100 moment_dev: addition of tests for complex commit 6c2630be4ce1a03fa583ae0487343c9ddde2bdef Author: Vandenplas, Jeremie Date: Fri Feb 28 23:33:00 2020 +0100 moment_dev: test progress commit 1743f87c83762b4bc053c6695756f9fb4464bfbd Author: Vandenplas, Jeremie Date: Fri Feb 28 22:25:43 2020 +0100 moment_dev: correction of moment for complex commit dfbdd9761560fa8891dbe18ffbac284986362e64 Merge: 195c62d 5f35833 Author: Vandenplas, Jeremie Date: Fri Feb 28 20:12:35 2020 +0100 Merge remote-tracking branch 'upstream/master' into moment_dev commit 195c62d7f16f8b25bda0bb6f9fce4565def74173 Author: Vandenplas, Jeremie Date: Sun Feb 23 20:26:09 2020 +0100 moment_dev: removed some explicit statements commit 7c5233f74ceb3d5a5637fba5f09923e49a056e4f Author: Vandenplas, Jeremie Date: Sat Feb 22 16:51:06 2020 +0100 moment_dev: addtition of test commit 5ff096836e32f44d297b0f8bb75f26a20902f5f6 Author: Vandenplas, Jeremie Date: Sat Feb 22 16:09:54 2020 +0100 moment_dev: addition of a function for central moment --- src/CMakeLists.txt | 1 + src/Makefile.manual | 6 + src/stdlib_experimental_stats.fypp | 103 +++- src/stdlib_experimental_stats.md | 86 ++- src/stdlib_experimental_stats_moment.fypp | 261 ++++++++ src/tests/stats/CMakeLists.txt | 1 + src/tests/stats/Makefile.manual | 2 +- src/tests/stats/test_moment.f90 | 706 ++++++++++++++++++++++ 8 files changed, 1158 insertions(+), 8 deletions(-) create mode 100644 src/stdlib_experimental_stats_moment.fypp create mode 100644 src/tests/stats/test_moment.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6b127c0b5..043e6932d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ set(fppFiles stdlib_experimental_optval.fypp stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.fypp + stdlib_experimental_stats_moment.fypp stdlib_experimental_stats_var.fypp stdlib_experimental_quadrature.fypp stdlib_experimental_quadrature_trapz.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index 9bffb4337..8e912741b 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -8,6 +8,7 @@ SRC = f18estop.f90 \ stdlib_experimental_quadrature_trapz.f90 \ stdlib_experimental_stats.f90 \ stdlib_experimental_stats_mean.f90 \ + stdlib_experimental_stats_moment.f90 \ stdlib_experimental_stats_var.f90 LIB = libstdlib.a @@ -46,6 +47,10 @@ stdlib_experimental_stats_mean.o: \ stdlib_experimental_optval.o \ stdlib_experimental_kinds.o \ stdlib_experimental_stats.o +stdlib_experimental_stats_moment.o: \ + stdlib_experimental_optval.o \ + stdlib_experimental_kinds.o \ + stdlib_experimental_stats.o stdlib_experimental_stats_var.o: \ stdlib_experimental_optval.o \ stdlib_experimental_kinds.o \ @@ -56,4 +61,5 @@ stdlib_experimental_io.f90: stdlib_experimental_io.fypp stdlib_experimental_quadrature.f90: stdlib_experimental_quadrature.fypp stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp +stdlib_experimental_stats_moment.f90: stdlib_experimental_stats_moment.fypp stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index c74f6ff6f..d4b971ae1 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -7,7 +7,7 @@ module stdlib_experimental_stats implicit none private ! Public API - public :: mean, var + public :: mean, moment, var interface mean #:for k1, t1 in RC_KINDS_TYPES @@ -209,5 +209,106 @@ module stdlib_experimental_stats end interface var + interface moment + #: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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in), optional :: mask + ${t1}$ :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment_all",rank, t1, k1, 'dp') + module function ${RName}$(x, order, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in), optional :: mask + real(dp) :: res + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in), optional :: mask + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment_mask_all",rank, t1, k1) + module function ${RName}$(x, order, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp') + module function ${RName}$(x, order, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment_mask",rank, t1, k1, 'dp') + module function ${RName}$(x, order, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + end interface moment end module stdlib_experimental_stats diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index 85b07b58f..debbd98a4 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -1,9 +1,29 @@ # Descriptive statistics -## Implemented - * `mean` - * `var` +## Implemented + + +* [`mean` - mean of array elements](#mean---mean-of-array-elements) + * [Description](#description) + * [Syntax](#syntax) + * [Arguments](#arguments) + * [Return value](#return-value) + * [Example](#example) +* [`moment` - central moment of array elements](#moment---central-moment-of-array-elements) + * [Description](#description-1) + * [Syntax](#syntax-1) + * [Arguments](#arguments-1) + * [Return value](#return-value-1) + * [Example](#example-1) +* [`var` - variance of array elements](#var---variance-of-array-elements) + * [Description](#description-2) + * [Syntax](#syntax-2) + * [Arguments](#arguments-2) + * [Return value](#return-value-2) + * [Example](#example-2) + + ## `mean` - mean of array elements @@ -28,7 +48,7 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a ### 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 `integer`, the result is of type `real(dp)`. If `dim` is absent, a scalar with the mean 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. @@ -49,6 +69,60 @@ program demo_mean end program demo_mean ``` +## `moment` - central moment 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`. + +The _k_-th order central moment is defined as : + +``` + moment(array) = 1/n sum_i (array(i) - mean(array))^k +``` + +where n is the number of elements. + +### Syntax + +`result = moment(array, order [, mask])` + +`result = moment(array, order, dim [, mask])` + +### Arguments + +`array`: Shall be an array of type `integer`, `real`, or `complex`. + +`order`: Shall be an scalar of type `integer`. + +`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to n, where n is the rank of `array`. + +`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`. + +### Return value + +If `array` is of type `real` 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 `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 + 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] +end program demo_moment +``` + ## `var` - variance of array elements ### Description @@ -58,7 +132,7 @@ Returns the variance of all the elements of `array`, or of the elements of `arra 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 + var(array) = 1/(n-1) sum_i (array(i) - mean(array))^2 ``` where n is the number of elements. @@ -108,7 +182,7 @@ program demo_var 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] + corrected=.false.) !returns [NaN, 0., 0.25] end program demo_var ``` diff --git a/src/stdlib_experimental_stats_moment.fypp b/src/stdlib_experimental_stats_moment.fypp new file mode 100644 index 000000000..c0ca9049d --- /dev/null +++ b/src/stdlib_experimental_stats_moment.fypp @@ -0,0 +1,261 @@ +#:include "common.fypp" +#:set RANKS = range(1, MAXRANK + 1) +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule (stdlib_experimental_stats) stdlib_experimental_stats_moment + + use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use stdlib_experimental_error, only: error_stop + use stdlib_experimental_optval, only: optval + implicit none + +contains + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("moment_all",rank, t1, k1) + module function ${RName}$(x, order, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + 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) + return + end if + + n = size(x, kind = int64) + mean = sum(x) / n + + res = sum((x - mean)**order) / n + + end function ${RName}$ + #:endfor + #:endfor + + + #: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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in), optional :: mask + real(dp) :: res + + real(dp) :: n, mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + n = size(x, kind = int64) + mean = sum(real(x, dp)) / n + + res = sum((real(x, dp) - mean)**order) / n + + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: 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')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + res = 0 + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = size(x, dim) + mean = sum(x, dim) / n + do i = 1, size(x, dim) + res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / n + + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(dp) :: n + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + res = 0 + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + 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)**order + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / n + + end function ${RName}$ + #:endfor + #:endfor + + + #: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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + 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 + + end function ${RName}$ + #:endfor + #:endfor + + + #: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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + + real(dp) :: n, mean + + n = count(mask, kind = int64) + mean = sum(real(x, dp), mask) / n + + res = sum((real(x, dp) - mean)**order, mask) / n + + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: 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')}$ + + res = 0 + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = count(mask, dim) + mean = sum(x, dim, mask) / n + 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 + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / n + + 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) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: order + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(dp) :: n${reduced_shape('x', rank, 'dim')}$ + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0 + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + 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)**order,& + 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / n + + end function ${RName}$ + #:endfor + #:endfor + +end submodule diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 8dc45370c..aa0a9e1ba 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,4 +1,5 @@ ADDTEST(mean) +ADDTEST(moment) ADDTEST(var) ADDTEST(varn) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 1d25a5509..aacaf98ab 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,3 @@ -PROGS_SRC = test_mean.f90 test_var.f90 +PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 include ../Makefile.manual.test.mk diff --git a/src/tests/stats/test_moment.f90 b/src/tests/stats/test_moment.f90 new file mode 100644 index 000000000..eaec6367d --- /dev/null +++ b/src/tests/stats/test_moment.f90 @@ -0,0 +1,706 @@ +program test_moment + use stdlib_experimental_error, only: assert + use stdlib_experimental_kinds, only: sp, dp, int32, int64 + use stdlib_experimental_stats, only: moment + 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_dp(d1, d) + call test_int32(int(d1, int32), int(d, int32)) + call test_int64(int(d1, int64), int(d, int64)) + + 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(:, :, :) + + order = 1 + + !1dim + print*,' test_sp_1dim', order + call assert( abs(moment(x1, order)) < sptol) + call assert( abs(moment(x1, order, dim=1)) < sptol) + + print*,' test_sp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_sp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5)) < sptol) + call assert( abs(moment(x1, order, 1, x1 < 5)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call assert( abs(moment(x2, order)) < sptol) + call assert( all( abs( moment(x2, order, 1)) < sptol)) + call assert( all( abs( moment(x2, order, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_sp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)) < sptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11)) < sptol)) + call assert( all( abs( moment(x2, order, 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; + + print*,' test_sp_3dim', order + call assert( abs(moment(x3, order)) < sptol) + call assert( all( abs( moment(x3, order, 1)) < sptol)) + call assert( all( abs( moment(x3, order, 2)) < sptol)) + call assert( all( abs( moment(x3, order, 3)) < sptol)) + + print*,' test_sp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_sp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) ) < sptol) + call assert( all( abs( moment(x3, order, 1, x3 < 45)) < sptol )) + call assert( all( abs( moment(x3, order, 2, x3 < 45)) < sptol )) + call assert( all( abs( moment(x3, order, 3, x3 < 45)) < sptol )) + + + order = 2 + + !1dim + print*,' test_sp_1dim', order + call assert( abs(moment(x1, order) - 2._sp) < sptol) + call assert( abs(moment(x1, order, dim=1) - 2._sp) < sptol) + + print*,' test_sp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_sp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5) - 1.25_sp) < sptol) + call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_sp) < sptol) + + !2dim + print*,' test_sp_2dim', order + call assert( abs(moment(x2, order) - 107.25_sp/9.) < sptol) + call assert( all( abs( moment(x2, order, 1) - [5._sp, 5._sp, 1.25_sp]) < sptol)) + call assert( all( abs( moment(x2, order, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]*2./3.) < sptol)) + + print*,' test_sp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_sp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)- 2.75_sp*3.) < sptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11) -& + [5._sp, 5._sp, 0.25_sp]) < sptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11) -& + [19._sp*2./3., 43._sp/9.*2., 0.25_sp , 0.25_sp]) < sptol)) + + !3dim + print*,' test_sp_3dim', order + call assert( abs(moment(x3, order) - 153.4_sp*35./36.) < sptol) + call assert( all( abs( moment(x3, order, 1) -& + reshape([20._sp / 3., 20._sp / 3., 5._sp / 3.,& + 4* 20._sp / 3., 4* 20._sp / 3., 4* 5._sp / 3.,& + 16* 20._sp / 3., 16* 20._sp / 3., 16* 5._sp / 3.],& + [size(x3,2), size(x3,3)])*3._sp/4.)& + < sptol)) + call assert( all( abs( moment(x3, order, 2) -& + reshape([19._sp, 43._sp / 3., 31._sp / 3. , 7.0_sp,& + 4* 19.0_sp, 4* 43._sp / 3., 4* 31._sp / 3. , 4* 7.0_sp,& + 16* 19.0_sp, 16* 43._sp / 3., 16* 31._sp / 3. , 16* 7.0_sp],& + [size(x3,1), size(x3,3)] )*2._sp/3.)& + < sptol)) + call assert( all( abs( moment(x3, order, 3) -& + reshape([ 7._sp/3., 21._sp, 175._sp/3.,& + 343._sp/3., 28._sp/3., 112._sp/3.,& + 84._sp, 448._sp/3., 189._sp,& + 700._sp/3., 847._sp/3., 336._sp],& + [size(x3,1), size(x3,2)] )*2./3.)& + < sptol)) + + print*,' test_sp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_sp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < sptol) + call assert( all( abs( moment(x3, order, 1, 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 assert( all( abs( moment(x3, order, 2, 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 assert( all( abs( moment(x3, order, 3, 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& + ], [size(x3,1), size(x3,2)] ))& + < sptol )) + + end subroutine + + subroutine test_dp(x1, x2) + real(dp), intent(in) :: x1(:), x2(:, :) + + integer :: order + real(dp), allocatable :: x3(:, :, :) + + order = 1 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order)) < dptol) + call assert( abs(moment(x1, order, dim=1)) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5)) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order)) < dptol) + call assert( all( abs( moment(x2, order, 1)) < dptol)) + call assert( all( abs( moment(x2, order, 2)) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + + !3dim + allocate(x3(size(x2,1),size(x2,2),3)) + x3(:,:,1)=x2; + x3(:,:,2)=x2*2; + x3(:,:,3)=x2*4; + + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order)) < dptol) + call assert( all( abs( moment(x3, order, 1)) < dptol)) + call assert( all( abs( moment(x3, order, 2)) < dptol)) + call assert( all( abs( moment(x3, order, 3)) < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) ) < dptol) + call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + + + order = 2 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order) - 2._dp) < dptol) + call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol) + call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2) -& + [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11) -& + [5._dp, 5._dp, 0.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11) -& + [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) + + !3dim + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol) + call assert( all( abs( moment(x3, order, 1) -& + reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,& + 4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,& + 16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],& + [size(x3,2), size(x3,3)])*3._dp/4.)& + < dptol)) + call assert( all( abs( moment(x3, order, 2) -& + reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],& + [size(x3,1), size(x3,3)] )*2._dp/3.)& + < dptol)) + call assert( all( abs( moment(x3, order, 3) -& + 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(x3,1), size(x3,2)] )*2./3.)& + < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) + call assert( all( abs( moment(x3, order, 1, 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 assert( all( abs( moment(x3, order, 2, 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 assert( all( abs( moment(x3, order, 3, 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& + ], [size(x3,1), size(x3,2)] ))& + < dptol )) + + end subroutine + + subroutine test_int32(x1, x2) + integer(int32), intent(in) :: x1(:), x2(:, :) + + integer :: order + integer(int32), allocatable :: x3(:, :, :) + + order = 1 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order)) < dptol) + call assert( abs(moment(x1, order, dim=1)) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5)) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order)) < dptol) + call assert( all( abs( moment(x2, order, 1)) < dptol)) + call assert( all( abs( moment(x2, order, 2)) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + + !3dim + allocate(x3(size(x2,1),size(x2,2),3)) + x3(:,:,1)=x2; + x3(:,:,2)=x2*2; + x3(:,:,3)=x2*4; + + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order)) < dptol) + call assert( all( abs( moment(x3, order, 1)) < dptol)) + call assert( all( abs( moment(x3, order, 2)) < dptol)) + call assert( all( abs( moment(x3, order, 3)) < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) ) < dptol) + call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + + + order = 2 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order) - 2._dp) < dptol) + call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol) + call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2) -& + [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11) -& + [5._dp, 5._dp, 0.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11) -& + [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) + + !3dim + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol) + call assert( all( abs( moment(x3, order, 1) -& + reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,& + 4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,& + 16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],& + [size(x3,2), size(x3,3)])*3._dp/4.)& + < dptol)) + call assert( all( abs( moment(x3, order, 2) -& + reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],& + [size(x3,1), size(x3,3)] )*2._dp/3.)& + < dptol)) + call assert( all( abs( moment(x3, order, 3) -& + 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(x3,1), size(x3,2)] )*2./3.)& + < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) + call assert( all( abs( moment(x3, order, 1, 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 assert( all( abs( moment(x3, order, 2, 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 assert( all( abs( moment(x3, order, 3, 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& + ], [size(x3,1), size(x3,2)] ))& + < dptol )) + + end subroutine + + subroutine test_int64(x1, x2) + integer(int64), intent(in) :: x1(:), x2(:, :) + + integer :: order + integer(int64), allocatable :: x3(:, :, :) + + order = 1 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order)) < dptol) + call assert( abs(moment(x1, order, dim=1)) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5)) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5)) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order)) < dptol) + call assert( all( abs( moment(x2, order, 1)) < dptol)) + call assert( all( abs( moment(x2, order, 2)) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11)) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11)) < dptol)) + + !3dim + allocate(x3(size(x2,1),size(x2,2),3)) + x3(:,:,1)=x2; + x3(:,:,2)=x2*2; + x3(:,:,3)=x2*4; + + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order)) < dptol) + call assert( all( abs( moment(x3, order, 1)) < dptol)) + call assert( all( abs( moment(x3, order, 2)) < dptol)) + call assert( all( abs( moment(x3, order, 3)) < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) ) < dptol) + call assert( all( abs( moment(x3, order, 1, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 2, x3 < 45)) < dptol )) + call assert( all( abs( moment(x3, order, 3, x3 < 45)) < dptol )) + + + order = 2 + + !1dim + print*,' test_dp_1dim', order + call assert( abs(moment(x1, order) - 2._dp) < dptol) + call assert( abs(moment(x1, order, dim=1) - 2._dp) < dptol) + + print*,' test_dp_1dim_mask', order + call assert( isnan(moment(x1, order, .false.))) + call assert( isnan(moment(x1, order, 1, .false.))) + + print*,' test_dp_1dim_mask_array', order + call assert( abs(moment(x1, order, x1 < 5) - 1.25_dp) < dptol) + call assert( abs(moment(x1, order, 1, x1 < 5) - 1.25_dp) < dptol) + + !2dim + print*,' test_dp_2dim', order + call assert( abs(moment(x2, order) - 107.25_dp/9.) < dptol) + call assert( all( abs( moment(x2, order, 1) - [5._dp, 5._dp, 1.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2) -& + [19._dp, 43._dp / 3., 31._dp / 3. , 7._dp]*2._dp/3.) < dptol)) + + print*,' test_dp_2dim_mask', order + call assert( isnan(moment(x2, order, .false.))) + call assert( any(isnan(moment(x2, order, 1, .false.)))) + call assert( any(isnan(moment(x2, order, 2, .false.)))) + + print*,' test_dp_2dim_mask_array', order + call assert( abs(moment(x2, order, x2 < 11)- 2.75_dp*3.) < dptol) + call assert( all( abs( moment(x2, order, 1, x2 < 11) -& + [5._dp, 5._dp, 0.25_dp]) < dptol)) + call assert( all( abs( moment(x2, order, 2, x2 < 11) -& + [19._dp*2./3., 43._dp/9.*2., 0.25_dp , 0.25_dp]) < dptol)) + + !3dim + print*,' test_dp_3dim', order + call assert( abs(moment(x3, order) - 153.4_dp*35./36.) < dptol) + call assert( all( abs( moment(x3, order, 1) -& + reshape([20._dp / 3., 20._dp / 3., 5._dp / 3.,& + 4* 20._dp / 3., 4* 20._dp / 3., 4* 5._dp / 3.,& + 16* 20._dp / 3., 16* 20._dp / 3., 16* 5._dp / 3.],& + [size(x3,2), size(x3,3)])*3._dp/4.)& + < dptol)) + call assert( all( abs( moment(x3, order, 2) -& + reshape([19._dp, 43._dp / 3., 31._dp / 3. , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3., 4* 31._dp / 3. , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3., 16* 31._dp / 3. , 16* 7.0_dp],& + [size(x3,1), size(x3,3)] )*2._dp/3.)& + < dptol)) + call assert( all( abs( moment(x3, order, 3) -& + 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(x3,1), size(x3,2)] )*2./3.)& + < dptol)) + + print*,' test_dp_3dim_mask', order + call assert( isnan(moment(x3, order, .false.))) + call assert( any(isnan(moment(x3, order, 1, .false.)))) + call assert( any(isnan(moment(x3, order, 2, .false.)))) + call assert( any(isnan(moment(x3, order, 3, .false.)))) + + print*,' test_dp_3dim_mask_array', order + call assert( abs(moment(x3, order, x3 < 11) - 7.7370242214532876_dp ) < dptol) + call assert( all( abs( moment(x3, order, 1, 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 assert( all( abs( moment(x3, order, 2, 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 assert( all( abs( moment(x3, order, 3, 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& + ], [size(x3,1), size(x3,2)] ))& + < dptol )) + + end subroutine + + subroutine test_csp(x1, x2) + complex(sp), intent(in) :: x1(:), x2(:, :) + + integer :: order + complex(sp), allocatable :: x3(:, :, :) + + order = 1 + + !1dim + print*,' test_sp_1dim', order + call assert( abs(moment(x1, order)) < sptol) + call assert( abs(moment(x1, order, dim=1)) < sptol) + + print*,' test_sp_1dim_mask', order + call assert( isnan(abs(moment(x1, order, .false.)))) + call assert( isnan(abs(moment(x1, order, 1, .false.)))) + + print*,' test_sp_1dim_mask_array', order + call assert( abs(moment(x1, order, aimag(x1) == 0)) < sptol) + call assert( abs(moment(x1, order, 1, aimag(x1) == 0)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call assert( abs(moment(x2, order)) < sptol) + call assert( all( abs( moment(x2, order, 1)) < sptol)) + call assert( all( abs( moment(x2, order, 2)) < sptol)) + + print*,' test_sp_2dim_mask', order + call assert( isnan(abs(moment(x2, order, .false.)))) + call assert( any(isnan(abs(moment(x2, order, 1, .false.))))) + call assert( any(isnan(abs(moment(x2, order, 2, .false.))))) + + print*,' test_sp_2dim_mask_array', order + call assert( abs(moment(x2, order, aimag(x2) == 0)) < sptol) + call assert( all( abs( moment(x2, order, 1, aimag(x2) == 0)) < sptol)) + call assert( any(isnan( abs( moment(x2, order, 2, aimag(x2) == 0))))) + + order = 2 + + !1dim + print*,' test_sp_1dim', order + call assert( abs(moment(x1, order) - (-6.459422410E-02,-0.556084037)) < sptol) + call assert( abs(moment(x1, order, dim=1) -& + (-6.459422410E-02,-0.556084037)) < sptol) + + print*,' test_sp_1dim_mask', order + call assert( isnan(abs(moment(x1, order, .false.)))) + call assert( isnan(abs(moment(x1, order, 1, .false.)))) + + print*,' test_sp_1dim_mask_array', order + call assert( abs(moment(x1, order, aimag(x1) == 0) -& + (8.969944715E-02,0.00000000)) < sptol) + call assert( abs(moment(x1, order, 1, aimag(x1) == 0) -& + (8.969944715E-02,0.00000000)) < sptol) + + !2dim + print*,' test_sp_2dim', order + call assert( abs(moment(x2, order) - (-0.163121477,-1.86906016)) < sptol) + call assert( all( abs( moment(x2, order, 1) -& + [(-6.459422410E-02,-0.556084037),& + (-0.581347823,-5.00475645),& + (-0.145336956,-1.25118911)]& + ) < sptol)) + call assert( all( abs( moment(x2, order, 2) -& + [(0.240498722,0.00000000),& + (-1.49895227,0.00000000),& + (1.15390968,0.00000000),& + (-0.569927275,0.00000000),& + (0.944405317,0.00000000)]& + ) < sptol)) + + print*,' test_sp_2dim_mask', order + call assert( isnan(abs(moment(x2, order, .false.)))) + call assert( any(isnan(abs(moment(x2, order, 1, .false.))))) + call assert( any(isnan(abs(moment(x2, order, 2, .false.))))) + + print*,' test_sp_2dim_mask_array', order + call assert( abs(moment(x2, order, aimag(x2) == 0)-& + (1.08109438,0.00000000)) < sptol) + call assert( all( abs( moment(x2, order, 1, aimag(x2)==0) -& + [(8.969944715E-02,0.00000000),& + (0.807295084,0.00000000),& + (0.201823771,0.00000000)]& + ) < sptol)) + + end subroutine +end program From 465a85d3106ffcb8b37eb0f8855069b338ca99f5 Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Mon, 9 Mar 2020 19:47:11 +0100 Subject: [PATCH 2/3] central_moment: n and mean computed outside select case --- src/stdlib_experimental_stats_moment.fypp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/stdlib_experimental_stats_moment.fypp b/src/stdlib_experimental_stats_moment.fypp index c0ca9049d..f36b3b5e3 100644 --- a/src/stdlib_experimental_stats_moment.fypp +++ b/src/stdlib_experimental_stats_moment.fypp @@ -82,12 +82,13 @@ contains return 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}$) - n = size(x, dim) - mean = sum(x, dim) / n do i = 1, size(x, dim) res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order end do @@ -121,12 +122,13 @@ contains return 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}$) - 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)**order end do @@ -198,12 +200,13 @@ contains real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$ ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + n = count(mask, dim) + mean = sum(x, dim, mask) / n + res = 0 select case(dim) #:for fi in range(1, rank+1) case(${fi}$) - n = count(mask, dim) - mean = sum(x, dim, mask) / n do i = 1, size(x, dim) res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**order,& #:if t1[0] == 'r' @@ -238,12 +241,13 @@ contains real(dp) :: n${reduced_shape('x', rank, 'dim')}$ real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + 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}$) - 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)**order,& 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) From 9ef4330165a85658aaa8fde3a48a33f633cb81cc Mon Sep 17 00:00:00 2001 From: "Vandenplas, Jeremie" Date: Thu, 12 Mar 2020 12:26:29 +0100 Subject: [PATCH 3/3] central_moment: changed mean to moment in the error messages --- src/stdlib_experimental_stats_moment.fypp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stdlib_experimental_stats_moment.fypp b/src/stdlib_experimental_stats_moment.fypp index f36b3b5e3..3a6d85dee 100644 --- a/src/stdlib_experimental_stats_moment.fypp +++ b/src/stdlib_experimental_stats_moment.fypp @@ -94,7 +94,7 @@ contains end do #:endfor case default - call error_stop("ERROR (mean): wrong dimension") + call error_stop("ERROR (moment): wrong dimension") end select res = res / n @@ -134,7 +134,7 @@ contains end do #:endfor case default - call error_stop("ERROR (mean): wrong dimension") + call error_stop("ERROR (moment): wrong dimension") end select res = res / n @@ -218,7 +218,7 @@ contains end do #:endfor case default - call error_stop("ERROR (mean): wrong dimension") + call error_stop("ERROR (moment): wrong dimension") end select res = res / n @@ -254,7 +254,7 @@ contains end do #:endfor case default - call error_stop("ERROR (mean): wrong dimension") + call error_stop("ERROR (moment): wrong dimension") end select res = res / n