Skip to content

Commit 2a0182a

Browse files
jvdp1jvdp1
committed
variance_dev: use fypp to avoid abs in real functions
1 parent da90a89 commit 2a0182a

File tree

2 files changed

+39
-208
lines changed

2 files changed

+39
-208
lines changed

‎src/stdlib_experimental_stats.fypp

Lines changed: 4 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -107,20 +107,9 @@ module stdlib_experimental_stats
107107

108108
interface var
109109

110-
#:for k1, t1 in REAL_KINDS_TYPES
110+
#:for k1, t1 in RC_KINDS_TYPES
111111
#:for rank in RANKS
112112
#:set RName = rname("var_all",rank, t1, k1)
113-
module function ${RName}$(x, mask) result(res)
114-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
115-
logical, intent(in), optional :: mask
116-
${t1}$ :: res
117-
end function ${RName}$
118-
#:endfor
119-
#:endfor
120-
121-
#:for k1, t1 in CMPLX_KINDS_TYPES
122-
#:for rank in RANKS
123-
#:set RName = rname("var_all",rank, t1, k1, 'r'+k1)
124113
module function ${RName}$(x, mask) result(res)
125114
${t1}$, intent(in) :: x${ranksuffix(rank)}$
126115
logical, intent(in), optional :: mask
@@ -140,21 +129,9 @@ module stdlib_experimental_stats
140129
#:endfor
141130
#:endfor
142131

143-
#:for k1, t1 in REAL_KINDS_TYPES
132+
#:for k1, t1 in RC_KINDS_TYPES
144133
#:for rank in RANKS
145134
#:set RName = rname("var",rank, t1, k1)
146-
module function ${RName}$(x, dim, mask) result(res)
147-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
148-
integer, intent(in) :: dim
149-
logical, intent(in), optional :: mask
150-
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
151-
end function ${RName}$
152-
#:endfor
153-
#:endfor
154-
155-
#:for k1, t1 in CMPLX_KINDS_TYPES
156-
#:for rank in RANKS
157-
#:set RName = rname("var",rank, t1, k1, 'r'+k1)
158135
module function ${RName}$(x, dim, mask) result(res)
159136
${t1}$, intent(in) :: x${ranksuffix(rank)}$
160137
integer, intent(in) :: dim
@@ -176,20 +153,9 @@ module stdlib_experimental_stats
176153
#:endfor
177154
#:endfor
178155

179-
#:for k1, t1 in REAL_KINDS_TYPES
156+
#:for k1, t1 in RC_KINDS_TYPES
180157
#:for rank in RANKS
181158
#:set RName = rname("var_mask_all",rank, t1, k1)
182-
module function ${RName}$(x, mask) result(res)
183-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
184-
logical, intent(in) :: mask${ranksuffix(rank)}$
185-
${t1}$ :: res
186-
end function ${RName}$
187-
#:endfor
188-
#:endfor
189-
190-
#:for k1, t1 in CMPLX_KINDS_TYPES
191-
#:for rank in RANKS
192-
#:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1)
193159
module function ${RName}$(x, mask) result(res)
194160
${t1}$, intent(in) :: x${ranksuffix(rank)}$
195161
logical, intent(in) :: mask${ranksuffix(rank)}$
@@ -209,21 +175,9 @@ module stdlib_experimental_stats
209175
#:endfor
210176
#:endfor
211177

212-
#:for k1, t1 in REAL_KINDS_TYPES
178+
#:for k1, t1 in RC_KINDS_TYPES
213179
#:for rank in RANKS
214180
#:set RName = rname("var_mask",rank, t1, k1)
215-
module function ${RName}$(x, dim, mask) result(res)
216-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
217-
integer, intent(in) :: dim
218-
logical, intent(in) :: mask${ranksuffix(rank)}$
219-
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
220-
end function ${RName}$
221-
#:endfor
222-
#:endfor
223-
224-
#:for k1, t1 in CMPLX_KINDS_TYPES
225-
#:for rank in RANKS
226-
#:set RName = rname("var_mask",rank, t1, k1, 'r'+k1)
227181
module function ${RName}$(x, dim, mask) result(res)
228182
${t1}$, intent(in) :: x${ranksuffix(rank)}$
229183
integer, intent(in) :: dim

‎src/stdlib_experimental_stats_var.fypp

Lines changed: 35 additions & 158 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#:include "common.fypp"
22
#:set RANKS = range(1, MAXRANK + 1)
3+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
34
submodule (stdlib_experimental_stats) stdlib_experimental_stats_var
45

56
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
@@ -9,55 +10,30 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_var
910

1011
contains
1112

12-
#:for k1, t1 in REAL_KINDS_TYPES
13+
#:for k1, t1 in RC_KINDS_TYPES
1314
#:for rank in RANKS
1415
#:set RName = rname("var_all",rank, t1, k1)
15-
module function ${RName}$(x, mask) result(res)
16-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
17-
logical, intent(in), optional :: mask
18-
${t1}$ :: res
19-
20-
${t1}$ :: n, mean
21-
22-
if (.not.optval(mask, .true.)) then
23-
res = ieee_value(1._${k1}$, ieee_quiet_nan)
24-
return
25-
end if
26-
27-
n = real(size(x, kind = int64), ${k1}$)
28-
mean = sum(x) / n
29-
30-
res = sum((x - mean)**2) / (n - 1._${k1}$)
31-
32-
end function ${RName}$
33-
#:endfor
34-
#:endfor
35-
36-
37-
#:for k1, t1 in CMPLX_KINDS_TYPES
38-
#:for rank in RANKS
39-
#:set RName = rname("var_all",rank, t1, k1, 'r'+k1)
4016
module function ${RName}$(x, mask) result(res)
4117
${t1}$, intent(in) :: x${ranksuffix(rank)}$
4218
logical, intent(in), optional :: mask
4319
real(${k1}$) :: res
4420

45-
real(${k1}$) :: n, mean
21+
real(${k1}$) :: n
22+
${t1}$ :: mean
4623

4724
if (.not.optval(mask, .true.)) then
4825
res = ieee_value(1._${k1}$, ieee_quiet_nan)
4926
return
5027
end if
5128

5229
n = real(size(x, kind = int64), ${k1}$)
30+
mean = sum(x) / n
5331

54-
!real part
55-
mean = sum(real(x)) / n
56-
res = sum((real(x) - mean)**2)
57-
58-
!imaginary part
59-
mean = sum(aimag(x)) / n
60-
res = (res + sum((aimag(x) - mean)**2)) / (n - 1._${k1}$)
32+
#:if t1[0] == 'r'
33+
res = sum((x - mean)**2) / (n - 1._${k1}$)
34+
#:else
35+
res = sum(abs(x - mean)**2) / (n - 1._${k1}$)
36+
#:endif
6137

6238
end function ${RName}$
6339
#:endfor
@@ -89,47 +65,9 @@ contains
8965
#:endfor
9066

9167

92-
#:for k1, t1 in REAL_KINDS_TYPES
68+
#:for k1, t1 in RC_KINDS_TYPES
9369
#:for rank in RANKS
9470
#:set RName = rname("var",rank, t1, k1)
95-
module function ${RName}$(x, dim, mask) result(res)
96-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
97-
integer, intent(in) :: dim
98-
logical, intent(in), optional :: mask
99-
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
100-
101-
integer :: i
102-
${t1}$ :: n
103-
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$
104-
105-
if (.not.optval(mask, .true.)) then
106-
res = ieee_value(1._${k1}$, ieee_quiet_nan)
107-
return
108-
end if
109-
110-
res = 0._${k1}$
111-
select case(dim)
112-
#:for fi in range(1, rank+1)
113-
case(${fi}$)
114-
n = real(size(x, dim), ${k1}$)
115-
mean = sum(x, dim) / n
116-
do i = 1, size(x, dim)
117-
res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2
118-
end do
119-
#:endfor
120-
case default
121-
call error_stop("ERROR (mean): wrong dimension")
122-
end select
123-
res = res / (n - 1._${k1}$)
124-
125-
end function ${RName}$
126-
#:endfor
127-
#:endfor
128-
129-
130-
#:for k1, t1 in CMPLX_KINDS_TYPES
131-
#:for rank in RANKS
132-
#:set RName = rname("var",rank, t1, k1, 'r'+k1)
13371
module function ${RName}$(x, dim, mask) result(res)
13472
${t1}$, intent(in) :: x${ranksuffix(rank)}$
13573
integer, intent(in) :: dim
@@ -138,7 +76,7 @@ contains
13876

13977
integer :: i
14078
real(${k1}$) :: n
141-
real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$
79+
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$
14280

14381
if (.not.optval(mask, .true.)) then
14482
res = ieee_value(1._${k1}$, ieee_quiet_nan)
@@ -150,15 +88,13 @@ contains
15088
#:for fi in range(1, rank+1)
15189
case(${fi}$)
15290
n = real(size(x, dim), ${k1}$)
153-
!real part
154-
mean = sum(real(x), dim) / n
155-
do i = 1, size(x, dim)
156-
res = res + (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2
157-
end do
158-
!imaginary part
159-
mean = sum(aimag(x), dim) / n
91+
mean = sum(x, dim) / n
16092
do i = 1, size(x, dim)
161-
res = res + (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2
93+
#:if t1[0] == 'r'
94+
res = res + (x${rankindice(':', 'i', rank, fi )}$ - mean)**2
95+
#:else
96+
res = res + abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2
97+
#:endif
16298
end do
16399
#:endfor
164100
case default
@@ -209,45 +145,25 @@ contains
209145
#:endfor
210146

211147

212-
#:for k1, t1 in REAL_KINDS_TYPES
148+
#:for k1, t1 in RC_KINDS_TYPES
213149
#:for rank in RANKS
214150
#:set RName = rname("var_mask_all",rank, t1, k1)
215-
module function ${RName}$(x, mask) result(res)
216-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
217-
logical, intent(in) :: mask${ranksuffix(rank)}$
218-
${t1}$ :: res
219-
220-
${t1}$ :: n, mean
221-
222-
n = real(count(mask, kind = int64), ${k1}$)
223-
mean = sum(x, mask) / n
224-
225-
res = sum((x - mean)**2, mask) / (n - 1._${k1}$)
226-
227-
end function ${RName}$
228-
#:endfor
229-
#:endfor
230-
231-
232-
#:for k1, t1 in CMPLX_KINDS_TYPES
233-
#:for rank in RANKS
234-
#:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1)
235151
module function ${RName}$(x, mask) result(res)
236152
${t1}$, intent(in) :: x${ranksuffix(rank)}$
237153
logical, intent(in) :: mask${ranksuffix(rank)}$
238154
real(${k1}$) :: res
239155

240-
real(${k1}$) :: n, mean
156+
real(${k1}$) :: n
157+
${t1}$ :: mean
241158

242159
n = real(count(mask, kind = int64), ${k1}$)
160+
mean = sum(x, mask) / n
243161

244-
!real part
245-
mean = sum(real(x), mask) / n
246-
res = sum((real(x) - mean)**2, mask)
247-
248-
!imaginary part
249-
mean = sum(aimag(x), mask) / n
250-
res = (res + sum((aimag(x) - mean)**2, mask)) / (n - 1._${k1}$)
162+
#:if t1[0] == 'r'
163+
res = sum((x - mean)**2, mask) / (n - 1._${k1}$)
164+
#:else
165+
res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$)
166+
#:endif
251167

252168
end function ${RName}$
253169
#:endfor
@@ -274,44 +190,9 @@ contains
274190
#:endfor
275191

276192

277-
#:for k1, t1 in REAL_KINDS_TYPES
193+
#:for k1, t1 in RC_KINDS_TYPES
278194
#:for rank in RANKS
279195
#:set RName = rname("var_mask",rank, t1, k1)
280-
module function ${RName}$(x, dim, mask) result(res)
281-
${t1}$, intent(in) :: x${ranksuffix(rank)}$
282-
integer, intent(in) :: dim
283-
logical, intent(in) :: mask${ranksuffix(rank)}$
284-
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
285-
286-
integer :: i
287-
${t1}$ :: n${reduced_shape('x', rank, 'dim')}$
288-
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$
289-
290-
res = 0._${k1}$
291-
select case(dim)
292-
#:for fi in range(1, rank+1)
293-
case(${fi}$)
294-
n = real(count(mask, dim), ${k1}$)
295-
mean = sum(x, dim, mask) / n
296-
do i = 1, size(x, dim)
297-
res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,&
298-
0._${k1}$,&
299-
mask${rankindice(':', 'i', rank, fi)}$)
300-
end do
301-
#:endfor
302-
case default
303-
call error_stop("ERROR (mean): wrong dimension")
304-
end select
305-
res = res / (n - 1._${k1}$)
306-
307-
end function ${RName}$
308-
#:endfor
309-
#:endfor
310-
311-
312-
#:for k1, t1 in CMPLX_KINDS_TYPES
313-
#:for rank in RANKS
314-
#:set RName = rname("var_mask",rank, t1, k1, 'r'+k1)
315196
module function ${RName}$(x, dim, mask) result(res)
316197
${t1}$, intent(in) :: x${ranksuffix(rank)}$
317198
integer, intent(in) :: dim
@@ -320,24 +201,20 @@ contains
320201

321202
integer :: i
322203
real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
323-
real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$
204+
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$
324205

325206
res = 0._${k1}$
326207
select case(dim)
327208
#:for fi in range(1, rank+1)
328209
case(${fi}$)
329210
n = real(count(mask, dim), ${k1}$)
330-
!real part
331-
mean = sum(real(x), dim, mask) / n
332-
do i = 1, size(x, dim)
333-
res = res + merge( (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,&
334-
0._${k1}$,&
335-
mask${rankindice(':', 'i', rank, fi)}$)
336-
end do
337-
!imaginary part
338-
mean = sum(aimag(x), dim, mask) / n
211+
mean = sum(x, dim, mask) / n
339212
do i = 1, size(x, dim)
340-
res = res + merge( (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,&
213+
#:if t1[0] == 'r'
214+
res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,&
215+
#:else
216+
res = res + merge( abs(x${rankindice(':', 'i', rank, fi )}$ - mean)**2,&
217+
#:endif
341218
0._${k1}$,&
342219
mask${rankindice(':', 'i', rank, fi)}$)
343220
end do

0 commit comments

Comments
 (0)