Skip to content

Commit 8dddff4

Browse files
jvdp1jvdp1
committed
progress
1 parent 4b96cdd commit 8dddff4

File tree

7 files changed

+396
-1
lines changed

7 files changed

+396
-1
lines changed

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ set(fppFiles
1818
stdlib_stats_corr.fypp
1919
stdlib_stats_cov.fypp
2020
stdlib_stats_mean.fypp
21+
stdlib_stats_median.fypp
2122
stdlib_stats_moment.fypp
2223
stdlib_stats_moment_all.fypp
2324
stdlib_stats_moment_mask.fypp

‎src/Makefile.manual

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ SRCFYPP =\
1818
stdlib_stats_corr.fypp \
1919
stdlib_stats_cov.fypp \
2020
stdlib_stats_mean.fypp \
21+
stdlib_stats_median.fypp \
2122
stdlib_stats_moment.fypp \
2223
stdlib_stats_moment_all.fypp \
2324
stdlib_stats_moment_mask.fypp \
@@ -108,6 +109,11 @@ stdlib_stats_mean.o: \
108109
stdlib_optval.o \
109110
stdlib_kinds.o \
110111
stdlib_stats.o
112+
stdlib_stats_median.o: \
113+
stdlib_optval.o \
114+
stdlib_kinds.o \
115+
stdlib_sorting.o \
116+
stdlib_stats.o
111117
stdlib_stats_moment.o: \
112118
stdlib_optval.o \
113119
stdlib_kinds.o \

‎src/common.fypp

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -147,4 +147,34 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
147147
#:endcall
148148
#:enddef
149149

150+
151+
#:def select_subvector(rank, idim)
152+
#:assert rank > 0
153+
#:call join_lines(joinstr=", ", prefix="(", suffix=")")
154+
#:for i in range(1, idim)
155+
${"j" + "_" * (i) }$
156+
#:endfor
157+
:
158+
#:for i in range(idim + 1, rank + 1)
159+
${"j" + "_" * (i) }$
160+
#:endfor
161+
#:endcall
162+
#:enddef
163+
164+
#:def reduce_subvector(rank, idim)
165+
#:assert rank > 0
166+
#:if rank > 1
167+
#:call join_lines(joinstr=", ", prefix="(", suffix=")")
168+
#:for i in range(1, idim)
169+
${"j" + "_" * (i) }$
170+
#:endfor
171+
#:for i in range(idim + 1, rank + 1)
172+
${"j" + "_" * (i) }$
173+
#:endfor
174+
#:endcall
175+
#:endif
176+
#:enddef
177+
178+
179+
150180
#:endmute

‎src/stdlib_stats.fypp

Lines changed: 57 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ module stdlib_stats
1111
implicit none
1212
private
1313
! Public API
14-
public :: corr, cov, mean, moment, var
14+
public :: corr, cov, mean, median, moment, var
1515

1616

1717
interface corr
@@ -304,6 +304,62 @@ module stdlib_stats
304304

305305
end interface mean
306306

307+
interface median
308+
!! version: experimental
309+
!!
310+
!! Median of array elements
311+
!! ([Specification](../page/specs/stdlib_stats.html#median))
312+
#:for k1, t1 in REAL_KINDS_TYPES
313+
#:for rank in RANKS
314+
#:set RName = rname("median_all",rank, t1, k1)
315+
module function ${RName}$ (x, mask) result(res)
316+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
317+
logical, intent(in), optional :: mask
318+
${t1}$ :: res
319+
end function ${RName}$
320+
#:endfor
321+
#:endfor
322+
323+
#:for k1, t1 in INT_KINDS_TYPES
324+
#:for rank in RANKS
325+
#:set RName = rname("median_all",rank, t1, k1,'dp')
326+
module function ${RName}$ (x, mask) result(res)
327+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
328+
logical, intent(in), optional :: mask
329+
real(dp) :: res
330+
end function ${RName}$
331+
#:endfor
332+
#:endfor
333+
334+
#:for k1, t1 in REAL_KINDS_TYPES
335+
#:for rank in RANKS
336+
#:set RName = rname("median",rank, t1, k1)
337+
module function ${RName}$(x, dim, mask) result(res)
338+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
339+
integer, intent(in) :: dim
340+
logical, intent(in), optional :: mask
341+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
342+
end function ${RName}$
343+
#:endfor
344+
#:endfor
345+
346+
#:for k1, t1 in INT_KINDS_TYPES
347+
#:for rank in RANKS
348+
#:set RName = rname("median",rank, t1, k1,'dp')
349+
module function ${RName}$(x, dim, mask) result(res)
350+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
351+
integer, intent(in) :: dim
352+
logical, intent(in), optional :: mask
353+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
354+
end function ${RName}$
355+
#:endfor
356+
#:endfor
357+
358+
359+
360+
361+
end interface
362+
307363

308364
interface var
309365
!! version: experimental

‎src/stdlib_stats_median.fypp

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
#:include "common.fypp"
2+
#:set RANKS = range(1, MAXRANK + 1)
3+
submodule (stdlib_stats) stdlib_stats_median
4+
5+
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
6+
use stdlib_error, only: error_stop
7+
use stdlib_optval, only: optval
8+
use stdlib_sorting, only: sort
9+
implicit none
10+
11+
contains
12+
13+
#:for k1, t1 in REAL_KINDS_TYPES
14+
#:for rank in RANKS
15+
#:set RName = rname("median_all",rank, t1, k1)
16+
module function ${RName}$ (x, mask) result(res)
17+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
18+
logical, intent(in), optional :: mask
19+
${t1}$ :: res
20+
21+
integer :: c, n
22+
${t1}$, allocatable :: x_tmp(:)
23+
24+
if (.not.optval(mask, .true.)) then
25+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
26+
return
27+
end if
28+
29+
x_tmp = reshape(x, [size(x)])
30+
31+
n = size(x_tmp)
32+
c = floor( (n + 1) / 2._${k1}$ )
33+
34+
call sort(x_tmp)
35+
36+
if (mod(n, 2) == 0) then
37+
res = sum(x_tmp(c:c+1)) / 2
38+
else
39+
res = x_tmp(c)
40+
endif
41+
42+
end function ${RName}$
43+
#:endfor
44+
#:endfor
45+
46+
#:for k1, t1 in INT_KINDS_TYPES
47+
#:for rank in RANKS
48+
#:set RName = rname("median_all",rank, t1, k1,'dp')
49+
module function ${RName}$ (x, mask) result(res)
50+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
51+
logical, intent(in), optional :: mask
52+
real(dp) :: res
53+
54+
integer :: c, n
55+
${t1}$, allocatable :: x_tmp(:)
56+
57+
if (.not.optval(mask, .true.)) then
58+
res = ieee_value(1._dp, ieee_quiet_nan)
59+
return
60+
end if
61+
62+
x_tmp = reshape(x, [size(x)])
63+
64+
n = size(x_tmp)
65+
c = floor( (n + 1) / 2.)
66+
67+
call sort(x_tmp)
68+
69+
if (mod(n, 2) == 0) then
70+
res = sum(real(x_tmp(c:c+1), kind=dp)) / 2._dp
71+
else
72+
res = x_tmp(c)
73+
endif
74+
75+
end function ${RName}$
76+
#:endfor
77+
#:endfor
78+
79+
#:for k1, t1 in REAL_KINDS_TYPES
80+
#:for rank in RANKS
81+
#:set RName = rname("median",rank, t1, k1)
82+
module function ${RName}$(x, dim, mask) result(res)
83+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
84+
integer, intent(in) :: dim
85+
logical, intent(in), optional :: mask
86+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
87+
88+
integer :: c, n
89+
#:for fj in range(1, rank+1)
90+
integer :: j${"_" * fj}$
91+
#:endfor
92+
${t1}$, allocatable :: x_tmp(:)
93+
94+
if (.not.optval(mask, .true.)) then
95+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
96+
return
97+
end if
98+
99+
n = size(x, dim)
100+
c = floor( (n + 1) / 2._${k1}$ )
101+
102+
allocate(x_tmp(n))
103+
104+
select case(dim)
105+
#:for fi in range(1, rank+1)
106+
case(${fi}$)
107+
#:for fj in range(1, fi)
108+
do j${"_" * fj}$ = 1, size(x, ${fj}$)
109+
#:endfor
110+
#:for fj in range(fi+1, rank+1)
111+
do j${"_" * fj}$ = 1, size(x, ${fj}$)
112+
#:endfor
113+
x_tmp(:) = x${select_subvector(rank, fi)}$
114+
call sort(x_tmp)
115+
116+
if (mod(n, 2) == 0) then
117+
res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._${k1}$
118+
else
119+
res${reduce_subvector(rank, fi)}$ = x_tmp(c)
120+
endif
121+
#:for fj in range(1, rank)
122+
end do
123+
#:endfor
124+
#:endfor
125+
case default
126+
call error_stop("ERROR (var): wrong dimension")
127+
end select
128+
129+
end function ${RName}$
130+
#:endfor
131+
#:endfor
132+
133+
#:for k1, t1 in INT_KINDS_TYPES
134+
#:for rank in RANKS
135+
#:set RName = rname("median",rank, t1, k1,'dp')
136+
module function ${RName}$(x, dim, mask) result(res)
137+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
138+
integer, intent(in) :: dim
139+
logical, intent(in), optional :: mask
140+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
141+
142+
integer :: c, n
143+
#:for fj in range(1, rank+1)
144+
integer :: j${"_" * fj}$
145+
#:endfor
146+
${t1}$, allocatable :: x_tmp(:)
147+
148+
if (.not.optval(mask, .true.)) then
149+
res = ieee_value(1._dp, ieee_quiet_nan)
150+
return
151+
end if
152+
153+
n = size(x, dim)
154+
c = floor( (n + 1) / 2._dp )
155+
156+
allocate(x_tmp(n))
157+
158+
select case(dim)
159+
#:for fi in range(1, rank+1)
160+
case(${fi}$)
161+
#:for fj in range(1, fi)
162+
do j${"_" * fj}$ = 1, size(x, ${fj}$)
163+
#:endfor
164+
#:for fj in range(fi+1, rank+1)
165+
do j${"_" * fj}$ = 1, size(x, ${fj}$)
166+
#:endfor
167+
x_tmp(:) = x${select_subvector(rank, fi)}$
168+
call sort(x_tmp)
169+
170+
if (mod(n, 2) == 0) then
171+
res${reduce_subvector(rank, fi)}$ = sum(x_tmp(c:c+1)) / 2._dp
172+
else
173+
res${reduce_subvector(rank, fi)}$ = x_tmp(c)
174+
endif
175+
#:for fj in range(1, rank)
176+
end do
177+
#:endfor
178+
#:endfor
179+
case default
180+
call error_stop("ERROR (var): wrong dimension")
181+
end select
182+
183+
end function ${RName}$
184+
#:endfor
185+
#:endfor
186+
187+
188+
end submodule

‎src/tests/stats/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
ADDTEST(corr)
22
ADDTEST(cov)
33
ADDTEST(mean)
4+
ADDTEST(median)
45
ADDTEST(moment)
56
ADDTEST(rawmoment)
67
ADDTEST(var)

0 commit comments

Comments
 (0)