Skip to content

Commit 71c30d1

Browse files
certikcertik
authored
Merge pull request #140 from fiolj/stats
Added support and test for mean of complex arrays
2 parents 2508db6 + 76bb81e commit 71c30d1

6 files changed

+141
-72
lines changed

‎src/common.fypp

+15
Original file line numberDiff line numberDiff line change
@@ -97,5 +97,20 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
9797
#:endif
9898
#:enddef
9999

100+
#! Generates a routine name from a generic name, rank, type and kind
101+
#!
102+
#! Args:
103+
#! gname (str): Generic name
104+
#! rank (integer): Rank if exist
105+
#! type (str): Type of the input
106+
#! kind (str): kind of inputs variable
107+
#! suffix (str): other identifier (could be used for output type/kind)
108+
#!
109+
#! Returns:
110+
#! A string with a new name
111+
#!
112+
#:def rname(gname, rank, type, kind, suffix='')
113+
$:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix)
114+
#:enddef
100115

101116
#:endmute

‎src/stdlib_experimental_stats.fypp

+29-28
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
#:include "common.fypp"
2-
32
#:set RANKS = range(1, MAXRANK + 1)
4-
5-
3+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
64
module stdlib_experimental_stats
75
use stdlib_experimental_kinds, only: sp, dp, qp, &
86
int8, int16, int32, int64
@@ -12,92 +10,95 @@ module stdlib_experimental_stats
1210
public :: mean
1311

1412
interface mean
15-
16-
#:for k1, t1 in REAL_KINDS_TYPES
13+
#:for k1, t1 in RC_KINDS_TYPES
1714
#:for rank in RANKS
18-
module function mean_${rank}$_all_${k1}$_${k1}$(x, mask) result(res)
15+
#:set RName = rname("mean_all",rank, t1, k1)
16+
module function ${RName}$ (x, mask) result(res)
1917
${t1}$, intent(in) :: x${ranksuffix(rank)}$
2018
logical, intent(in), optional :: mask
2119
${t1}$ :: res
22-
end function mean_${rank}$_all_${k1}$_${k1}$
20+
end function ${RName}$
2321
#:endfor
2422
#:endfor
2523

2624
#:for k1, t1 in INT_KINDS_TYPES
2725
#:for rank in RANKS
28-
module function mean_${rank}$_all_${k1}$_dp(x, mask) result(res)
26+
#:set RName = rname('mean_all', rank, t1, k1,'dp')
27+
module function ${RName}$(x, mask) result(res)
2928
${t1}$, intent(in) :: x${ranksuffix(rank)}$
3029
logical, intent(in), optional :: mask
3130
real(dp) :: res
32-
end function mean_${rank}$_all_${k1}$_dp
31+
end function ${RName}$
3332
#:endfor
3433
#:endfor
3534

36-
#:for k1, t1 in REAL_KINDS_TYPES
35+
#:for k1, t1 in RC_KINDS_TYPES
3736
#:for rank in RANKS
38-
module function mean_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res)
37+
#:set RName = rname("mean",rank, t1, k1)
38+
module function ${RName}$(x, dim, mask) result(res)
3939
${t1}$, intent(in) :: x${ranksuffix(rank)}$
4040
integer, intent(in) :: dim
4141
logical, intent(in), optional :: mask
4242
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
43-
end function mean_${rank}$_${k1}$_${k1}$
43+
end function ${RName}$
4444
#:endfor
4545
#:endfor
4646

4747
#:for k1, t1 in INT_KINDS_TYPES
4848
#:for rank in RANKS
49-
module function mean_${rank}$_${k1}$_dp(x, dim, mask) result(res)
49+
#:set RName = rname("mean",rank, t1, k1,'dp')
50+
module function ${RName}$(x, dim, mask) result(res)
5051
${t1}$, intent(in) :: x${ranksuffix(rank)}$
5152
integer, intent(in) :: dim
5253
logical, intent(in), optional :: mask
5354
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
54-
end function mean_${rank}$_${k1}$_dp
55+
end function ${RName}$
5556
#:endfor
5657
#:endfor
5758

58-
59-
#:for k1, t1 in REAL_KINDS_TYPES
59+
#:for k1, t1 in RC_KINDS_TYPES
6060
#:for rank in RANKS
61-
module function mean_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res)
61+
#:set RName = rname('mean_mask_all',rank, t1, k1)
62+
module function ${RName}$(x, mask) result(res)
6263
${t1}$, intent(in) :: x${ranksuffix(rank)}$
6364
logical, intent(in) :: mask${ranksuffix(rank)}$
6465
${t1}$ :: res
65-
end function mean_${rank}$_mask_all_${k1}$_${k1}$
66+
end function ${RName}$
6667
#:endfor
6768
#:endfor
6869

69-
7070
#:for k1, t1 in INT_KINDS_TYPES
7171
#:for rank in RANKS
72-
module function mean_${rank}$_mask_all_${k1}$_dp(x, mask) result(res)
72+
#:set RName = rname('mean_mask_all',rank, t1, k1, 'dp')
73+
module function ${RName}$(x, mask) result(res)
7374
${t1}$, intent(in) :: x${ranksuffix(rank)}$
7475
logical, intent(in) :: mask${ranksuffix(rank)}$
7576
real(dp) :: res
76-
end function mean_${rank}$_mask_all_${k1}$_dp
77+
end function ${RName}$
7778
#:endfor
7879
#:endfor
7980

80-
81-
#:for k1, t1 in REAL_KINDS_TYPES
81+
#:for k1, t1 in RC_KINDS_TYPES
8282
#:for rank in RANKS
83-
module function mean_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res)
83+
#:set RName = rname('mean_mask',rank, t1, k1)
84+
module function ${RName}$(x, dim, mask) result(res)
8485
${t1}$, intent(in) :: x${ranksuffix(rank)}$
8586
integer, intent(in) :: dim
8687
logical, intent(in) :: mask${ranksuffix(rank)}$
8788
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
88-
end function mean_${rank}$_mask_${k1}$_${k1}$
89+
end function ${RName}$
8990
#:endfor
9091
#:endfor
9192

92-
9393
#:for k1, t1 in INT_KINDS_TYPES
9494
#:for rank in RANKS
95-
module function mean_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res)
95+
#:set RName = rname('mean_mask',rank, t1, k1, 'dp')
96+
module function ${RName}$(x, dim, mask) result(res)
9697
${t1}$, intent(in) :: x${ranksuffix(rank)}$
9798
integer, intent(in) :: dim
9899
logical, intent(in) :: mask${ranksuffix(rank)}$
99100
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
100-
end function mean_${rank}$_mask_${k1}$_dp
101+
end function ${RName}$
101102
#:endfor
102103
#:endfor
103104

‎src/stdlib_experimental_stats.md

+2-2
Original file line numberDiff line numberDiff line change
@@ -18,15 +18,15 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
1818

1919
### Arguments
2020

21-
`array`: Shall be an array of type `integer`, or `real`.
21+
`array`: Shall be an array of type `integer`, `real`, or `complex`.
2222

2323
`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`.
2424

2525
`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`.
2626

2727
### Return value
2828

29-
If `array` is of type `real`, the result is of the same type as `array`.
29+
If `array` is of type `real` or `complex`, the result is of the same type as `array`.
3030
If `array` is of type `integer`, the result is of type `double precision`.
3131

3232
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.

‎src/stdlib_experimental_stats_mean.fypp

+31-28
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
#:include "common.fypp"
2-
32
#:set RANKS = range(1, MAXRANK + 1)
4-
5-
3+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
64
submodule (stdlib_experimental_stats) stdlib_experimental_stats_mean
75

86
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
@@ -12,28 +10,29 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_mean
1210

1311
contains
1412

15-
#:for k1, t1 in REAL_KINDS_TYPES
13+
#:for k1, t1 in RC_KINDS_TYPES
1614
#:for rank in RANKS
17-
module function mean_${rank}$_all_${k1}$_${k1}$(x, mask) result(res)
15+
#:set RName = rname("mean_all",rank, t1, k1)
16+
module function ${RName}$ (x, mask) result(res)
1817
${t1}$, intent(in) :: x${ranksuffix(rank)}$
1918
logical, intent(in), optional :: mask
2019
${t1}$ :: res
2120

2221
if (.not.optval(mask, .true.)) then
23-
res = ieee_value(res, ieee_quiet_nan)
22+
res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan)
2423
return
2524
end if
2625

2726
res = sum(x) / real(size(x, kind = int64), ${k1}$)
2827

29-
end function mean_${rank}$_all_${k1}$_${k1}$
28+
end function ${RName}$
3029
#:endfor
3130
#:endfor
3231

33-
3432
#:for k1, t1 in INT_KINDS_TYPES
3533
#:for rank in RANKS
36-
module function mean_${rank}$_all_${k1}$_dp(x, mask) result(res)
34+
#:set RName = rname('mean_all', rank, t1, k1,'dp')
35+
module function ${RName}$(x, mask) result(res)
3736
${t1}$, intent(in) :: x${ranksuffix(rank)}$
3837
logical, intent(in), optional :: mask
3938
real(dp) :: res
@@ -45,21 +44,22 @@ contains
4544

4645
res = sum(real(x, dp)) / real(size(x, kind = int64), dp)
4746

48-
end function mean_${rank}$_all_${k1}$_dp
47+
end function ${RName}$
4948
#:endfor
5049
#:endfor
5150

5251

53-
#:for k1, t1 in REAL_KINDS_TYPES
52+
#:for k1, t1 in RC_KINDS_TYPES
5453
#:for rank in RANKS
55-
module function mean_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res)
54+
#:set RName = rname("mean",rank, t1, k1)
55+
module function ${RName}$(x, dim, mask) result(res)
5656
${t1}$, intent(in) :: x${ranksuffix(rank)}$
5757
integer, intent(in) :: dim
5858
logical, intent(in), optional :: mask
5959
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
6060

6161
if (.not.optval(mask, .true.)) then
62-
res = ieee_value(res, ieee_quiet_nan)
62+
res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan)
6363
return
6464
end if
6565

@@ -69,14 +69,15 @@ contains
6969
call error_stop("ERROR (mean): wrong dimension")
7070
end if
7171

72-
end function mean_${rank}$_${k1}$_${k1}$
72+
end function ${RName}$
7373
#:endfor
7474
#:endfor
7575

7676

7777
#:for k1, t1 in INT_KINDS_TYPES
7878
#:for rank in RANKS
79-
module function mean_${rank}$_${k1}$_dp(x, dim, mask) result(res)
79+
#:set RName = rname("mean",rank, t1, k1,'dp')
80+
module function ${RName}$(x, dim, mask) result(res)
8081
${t1}$, intent(in) :: x${ranksuffix(rank)}$
8182
integer, intent(in) :: dim
8283
logical, intent(in), optional :: mask
@@ -93,42 +94,43 @@ contains
9394
call error_stop("ERROR (mean): wrong dimension")
9495
end if
9596

96-
end function mean_${rank}$_${k1}$_dp
97+
end function ${RName}$
9798
#:endfor
9899
#:endfor
99100

100-
101-
#:for k1, t1 in REAL_KINDS_TYPES
101+
#:for k1, t1 in RC_KINDS_TYPES
102102
#:for rank in RANKS
103-
module function mean_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res)
103+
#:set RName = rname('mean_mask_all',rank, t1, k1)
104+
module function ${RName}$(x, mask) result(res)
104105
${t1}$, intent(in) :: x${ranksuffix(rank)}$
105106
logical, intent(in) :: mask${ranksuffix(rank)}$
106107
${t1}$ :: res
107108

108109
res = sum(x, mask) / real(count(mask, kind = int64), ${k1}$)
109110

110-
end function mean_${rank}$_mask_all_${k1}$_${k1}$
111+
end function ${RName}$
111112
#:endfor
112113
#:endfor
113114

114115

115116
#:for k1, t1 in INT_KINDS_TYPES
116117
#:for rank in RANKS
117-
module function mean_${rank}$_mask_all_${k1}$_dp(x, mask) result(res)
118+
#:set RName = rname('mean_mask_all',rank, t1, k1, 'dp')
119+
module function ${RName}$(x, mask) result(res)
118120
${t1}$, intent(in) :: x${ranksuffix(rank)}$
119121
logical, intent(in) :: mask${ranksuffix(rank)}$
120122
real(dp) :: res
121123

122124
res = sum(real(x, dp), mask) / real(count(mask, kind = int64), dp)
123125

124-
end function mean_${rank}$_mask_all_${k1}$_dp
126+
end function ${RName}$
125127
#:endfor
126128
#:endfor
127129

128-
129-
#:for k1, t1 in REAL_KINDS_TYPES
130+
#:for k1, t1 in RC_KINDS_TYPES
130131
#:for rank in RANKS
131-
module function mean_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res)
132+
#:set RName = rname('mean_mask',rank, t1, k1)
133+
module function ${RName}$(x, dim, mask) result(res)
132134
${t1}$, intent(in) :: x${ranksuffix(rank)}$
133135
integer, intent(in) :: dim
134136
logical, intent(in) :: mask${ranksuffix(rank)}$
@@ -140,14 +142,15 @@ contains
140142
call error_stop("ERROR (mean): wrong dimension")
141143
end if
142144

143-
end function mean_${rank}$_mask_${k1}$_${k1}$
145+
end function ${RName}$
144146
#:endfor
145147
#:endfor
146148

147149

148150
#:for k1, t1 in INT_KINDS_TYPES
149151
#:for rank in RANKS
150-
module function mean_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res)
152+
#:set RName = rname('mean_mask',rank, t1, k1, 'dp')
153+
module function ${RName}$(x, dim, mask) result(res)
151154
${t1}$, intent(in) :: x${ranksuffix(rank)}$
152155
integer, intent(in) :: dim
153156
logical, intent(in) :: mask${ranksuffix(rank)}$
@@ -159,7 +162,7 @@ contains
159162
call error_stop("ERROR (mean): wrong dimension")
160163
end if
161164

162-
end function mean_${rank}$_mask_${k1}$_dp
165+
end function ${RName}$
163166
#:endfor
164167
#:endfor
165168

‎src/tests/stats/test_mean.f90

+33
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,13 @@ program test_mean
1313
real(sp), allocatable :: s(:, :)
1414
real(dp), allocatable :: d(:, :)
1515

16+
complex(dp), allocatable :: cs(:, :)
17+
complex(dp), allocatable :: cd(:, :)
18+
1619
real(dp), allocatable :: d3(:, :, :)
1720
real(dp), allocatable :: d4(:, :, :, :)
1821

22+
complex(dp), allocatable :: cd3(:, :, :)
1923

2024
!sp
2125
call loadtxt("array3.dat", s)
@@ -36,6 +40,23 @@ program test_mean
3640
call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) < dptol)
3741
call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) < dptol)
3842

43+
!csp
44+
45+
call loadtxt("array3.dat", d)
46+
cs = cmplx(1._sp, 1._sp)*d
47+
48+
call assert( abs(mean(cs) - sum(cs)/real(size(cs), sp)) < sptol)
49+
call assert( sum( abs( mean(cs,1) - sum(cs,1)/real(size(cs,1), sp) )) < sptol)
50+
call assert( sum( abs( mean(cs,2) - sum(cs,2)/real(size(cs,2), sp) )) < sptol)
51+
52+
!cdp
53+
54+
call loadtxt("array3.dat", d)
55+
cd = cmplx(1._dp, 1._dp)*d
56+
57+
call assert( abs(mean(cd) - sum(cd)/real(size(cd), dp)) < dptol)
58+
call assert( sum( abs( mean(cd,1) - sum(cd,1)/real(size(cd,1), dp) )) < dptol)
59+
call assert( sum( abs( mean(cd,2) - sum(cd,2)/real(size(cd,2), dp) )) < dptol)
3960

4061
! check mask = .false.
4162

@@ -75,6 +96,18 @@ program test_mean
7596
call assert( sum( abs( mean(d3,2) - sum(d3,2)/real(size(d3,2), dp) )) < dptol)
7697
call assert( sum( abs( mean(d3,3) - sum(d3,3)/real(size(d3,3), dp) )) < dptol)
7798

99+
!cdp rank 3
100+
allocate(cd3(size(d,1),size(d,2),3))
101+
cd3(:,:,1)=d;
102+
cd3(:,:,2)=d*1.5;
103+
cd3(:,:,3)=d*4;
104+
cd3 = cmplx(1._sp, 1._sp)*cd3
105+
106+
call assert( abs(mean(cd3) - sum(cd3)/real(size(cd3), dp)) < dptol)
107+
call assert( sum( abs( mean(cd3,1) - sum(cd3,1)/real(size(cd3,1), dp) )) < dptol)
108+
call assert( sum( abs( mean(cd3,2) - sum(cd3,2)/real(size(cd3,2), dp) )) < dptol)
109+
call assert( sum( abs( mean(cd3,3) - sum(cd3,3)/real(size(cd3,3), dp) )) < dptol)
110+
78111

79112
!dp rank 4
80113
allocate(d4(size(d,1),size(d,2),3,9))

0 commit comments

Comments
 (0)