Skip to content

Commit 95a7771

Browse files
committed
mv_to_central_moment: simplification
1 parent c37124a commit 95a7771

File tree

1 file changed

+20
-78
lines changed

1 file changed

+20
-78
lines changed

src/stdlib_experimental_stats_moment.fypp

Lines changed: 20 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -82,28 +82,16 @@ contains
8282
logical, intent(in), optional :: mask
8383
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
8484

85-
integer :: i
86-
real(${k1}$) :: n
87-
8885
if (.not.optval(mask, .true.)) then
8986
res = ieee_value(1._${k1}$, ieee_quiet_nan)
9087
return
9188
end if
9289

93-
n = size(x, dim)
94-
95-
res = 0
96-
select case(dim)
97-
#:for fi in range(1, rank+1)
98-
case(${fi}$)
99-
do i = 1, size(x, ${fi}$)
100-
res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - center)**order
101-
end do
102-
#:endfor
103-
case default
104-
call error_stop("ERROR (moment): wrong dimension")
105-
end select
106-
res = res / n
90+
if (dim >= 1 .and. dim <= ${rank}$) then
91+
res = sum((x - center)**order, dim) / size(x, dim)
92+
else
93+
call error_stop("ERROR (moment): wrong dimension")
94+
end if
10795

10896
end function ${RName}$
10997
#:endfor
@@ -169,29 +157,16 @@ contains
169157
logical, intent(in), optional :: mask
170158
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
171159

172-
integer :: i
173-
real(dp) :: n
174-
175160
if (.not.optval(mask, .true.)) then
176161
res = ieee_value(1._dp, ieee_quiet_nan)
177162
return
178163
end if
179164

180-
n = size(x, dim)
181-
182-
res = 0
183-
select case(dim)
184-
#:for fi in range(1, rank+1)
185-
case(${fi}$)
186-
do i = 1, size(x, ${fi}$)
187-
res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -&
188-
center)**order
189-
end do
190-
#:endfor
191-
case default
192-
call error_stop("ERROR (moment): wrong dimension")
193-
end select
194-
res = res / n
165+
if (dim >= 1 .and. dim <= ${rank}$) then
166+
res = sum( (real(x, dp) - center)**order, dim) / size(x, dim)
167+
else
168+
call error_stop("ERROR (moment): wrong dimension")
169+
end if
195170

196171
end function ${RName}$
197172
#:endfor
@@ -308,30 +283,11 @@ contains
308283
logical, intent(in) :: mask${ranksuffix(rank)}$
309284
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
310285

311-
integer :: i
312-
real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
313-
314-
n = count(mask, dim)
315-
316-
res = 0
317-
select case(dim)
318-
#:for fi in range(1, rank+1)
319-
case(${fi}$)
320-
do i = 1, size(x, ${fi}$)
321-
res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ -&
322-
center)**order,&
323-
#:if t1[0] == 'r'
324-
0._${k1}$,&
325-
#:else
326-
cmplx(0,0,kind=${k1}$),&
327-
#:endif
328-
mask${select_subarray(rank, [(fi, 'i')])}$)
329-
end do
330-
#:endfor
331-
case default
332-
call error_stop("ERROR (moment): wrong dimension")
333-
end select
334-
res = res / n
286+
if (dim >= 1 .and. dim <= ${rank}$) then
287+
res = sum((x - center)**order, dim, mask) / count(mask, dim)
288+
else
289+
call error_stop("ERROR (moment): wrong dimension")
290+
end if
335291

336292
end function ${RName}$
337293
#:endfor
@@ -405,25 +361,11 @@ contains
405361
logical, intent(in) :: mask${ranksuffix(rank)}$
406362
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
407363

408-
integer :: i
409-
real(dp) :: n${reduced_shape('x', rank, 'dim')}$
410-
411-
n = count(mask, dim)
412-
413-
res = 0
414-
select case(dim)
415-
#:for fi in range(1, rank+1)
416-
case(${fi}$)
417-
do i = 1, size(x, ${fi}$)
418-
res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -&
419-
center)**order,&
420-
0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
421-
end do
422-
#:endfor
423-
case default
424-
call error_stop("ERROR (moment): wrong dimension")
425-
end select
426-
res = res / n
364+
if (dim >= 1 .and. dim <= ${rank}$) then
365+
res = sum(( real(x, dp) - center)**order, dim, mask) / count(mask, dim)
366+
else
367+
call error_stop("ERROR (moment): wrong dimension")
368+
end if
427369

428370
end function ${RName}$
429371
#:endfor

0 commit comments

Comments
 (0)