Skip to content

Commit c37124a

Browse files
jvdp1jvdp1
committed
mv_to_central_moment: addition of a function with scalar center
1 parent cd2c16c commit c37124a

File tree

3 files changed

+349
-8
lines changed

3 files changed

+349
-8
lines changed

‎src/stdlib_experimental_stats.fypp

Lines changed: 57 additions & 0 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 REDRANKS = range(2, MAXRANK + 1)
34
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
45
module stdlib_experimental_stats
56
use stdlib_experimental_kinds, only: sp, dp, qp, &
@@ -236,6 +237,20 @@ module stdlib_experimental_stats
236237
#:endfor
237238
#:endfor
238239

240+
#:for k1, t1 in RC_KINDS_TYPES
241+
#:for rank in REDRANKS
242+
#:set RName = rname("moment_scalar",rank, t1, k1)
243+
module function ${RName}$(x, order, dim, center, mask) result(res)
244+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
245+
integer, intent(in) :: order
246+
integer, intent(in) :: dim
247+
${t1}$, intent(in) :: center
248+
logical, intent(in), optional :: mask
249+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
250+
end function ${RName}$
251+
#:endfor
252+
#:endfor
253+
239254
#:for k1, t1 in RC_KINDS_TYPES
240255
#:for rank in RANKS
241256
#:set RName = rname("moment",rank, t1, k1)
@@ -250,6 +265,20 @@ module stdlib_experimental_stats
250265
#:endfor
251266
#:endfor
252267

268+
#:for k1, t1 in INT_KINDS_TYPES
269+
#:for rank in REDRANKS
270+
#:set RName = rname("moment_scalar",rank, t1, k1, 'dp')
271+
module function ${RName}$(x, order, dim, center, mask) result(res)
272+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
273+
integer, intent(in) :: order
274+
integer, intent(in) :: dim
275+
real(dp),intent(in) :: center
276+
logical, intent(in), optional :: mask
277+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
278+
end function ${RName}$
279+
#:endfor
280+
#:endfor
281+
253282
#:for k1, t1 in INT_KINDS_TYPES
254283
#:for rank in RANKS
255284
#:set RName = rname("moment",rank, t1, k1, 'dp')
@@ -290,6 +319,20 @@ module stdlib_experimental_stats
290319
#:endfor
291320
#:endfor
292321

322+
#:for k1, t1 in RC_KINDS_TYPES
323+
#:for rank in REDRANKS
324+
#:set RName = rname("moment_mask_scalar",rank, t1, k1)
325+
module function ${RName}$(x, order, dim, center, mask) result(res)
326+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
327+
integer, intent(in) :: order
328+
integer, intent(in) :: dim
329+
${t1}$, intent(in) :: center
330+
logical, intent(in) :: mask${ranksuffix(rank)}$
331+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
332+
end function ${RName}$
333+
#:endfor
334+
#:endfor
335+
293336
#:for k1, t1 in RC_KINDS_TYPES
294337
#:for rank in RANKS
295338
#:set RName = rname("moment_mask",rank, t1, k1)
@@ -304,6 +347,20 @@ module stdlib_experimental_stats
304347
#:endfor
305348
#:endfor
306349

350+
#:for k1, t1 in INT_KINDS_TYPES
351+
#:for rank in REDRANKS
352+
#:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp')
353+
module function ${RName}$(x, order, dim, center, mask) result(res)
354+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
355+
integer, intent(in) :: order
356+
integer, intent(in) :: dim
357+
real(dp), intent(in) :: center
358+
logical, intent(in) :: mask${ranksuffix(rank)}$
359+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
360+
end function ${RName}$
361+
#:endfor
362+
#:endfor
363+
307364
#:for k1, t1 in INT_KINDS_TYPES
308365
#:for rank in RANKS
309366
#:set RName = rname("moment_mask",rank, t1, k1, 'dp')

‎src/stdlib_experimental_stats_moment.fypp

Lines changed: 157 additions & 0 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 REDRANKS = range(2, MAXRANK + 1)
34
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
45
submodule (stdlib_experimental_stats) stdlib_experimental_stats_moment
56

@@ -70,6 +71,45 @@ contains
7071
#:endfor
7172

7273

74+
#:for k1, t1 in RC_KINDS_TYPES
75+
#:for rank in REDRANKS
76+
#:set RName = rname("moment_scalar",rank, t1, k1)
77+
module function ${RName}$(x, order, dim, center, mask) result(res)
78+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
79+
integer, intent(in) :: order
80+
integer, intent(in) :: dim
81+
${t1}$, intent(in) :: center
82+
logical, intent(in), optional :: mask
83+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
84+
85+
integer :: i
86+
real(${k1}$) :: n
87+
88+
if (.not.optval(mask, .true.)) then
89+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
90+
return
91+
end if
92+
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
107+
108+
end function ${RName}$
109+
#:endfor
110+
#:endfor
111+
112+
73113
#:for k1, t1 in RC_KINDS_TYPES
74114
#:for rank in RANKS
75115
#:set RName = rname("moment",rank, t1, k1)
@@ -118,6 +158,46 @@ contains
118158
#:endfor
119159

120160

161+
#:for k1, t1 in INT_KINDS_TYPES
162+
#:for rank in REDRANKS
163+
#:set RName = rname("moment_scalar",rank, t1, k1, 'dp')
164+
module function ${RName}$(x, order, dim, center, mask) result(res)
165+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
166+
integer, intent(in) :: order
167+
integer, intent(in) :: dim
168+
real(dp),intent(in) :: center
169+
logical, intent(in), optional :: mask
170+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
171+
172+
integer :: i
173+
real(dp) :: n
174+
175+
if (.not.optval(mask, .true.)) then
176+
res = ieee_value(1._dp, ieee_quiet_nan)
177+
return
178+
end if
179+
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
195+
196+
end function ${RName}$
197+
#:endfor
198+
#:endfor
199+
200+
121201
#:for k1, t1 in INT_KINDS_TYPES
122202
#:for rank in RANKS
123203
#:set RName = rname("moment",rank, t1, k1, 'dp')
@@ -217,6 +297,47 @@ contains
217297
#:endfor
218298

219299

300+
#:for k1, t1 in RC_KINDS_TYPES
301+
#:for rank in REDRANKS
302+
#:set RName = rname("moment_mask_scalar",rank, t1, k1)
303+
module function ${RName}$(x, order, dim, center, mask) result(res)
304+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
305+
integer, intent(in) :: order
306+
integer, intent(in) :: dim
307+
${t1}$, intent(in) :: center
308+
logical, intent(in) :: mask${ranksuffix(rank)}$
309+
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
310+
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
335+
336+
end function ${RName}$
337+
#:endfor
338+
#:endfor
339+
340+
220341
#:for k1, t1 in RC_KINDS_TYPES
221342
#:for rank in RANKS
222343
#:set RName = rname("moment_mask",rank, t1, k1)
@@ -273,6 +394,42 @@ contains
273394
#:endfor
274395

275396

397+
#:for k1, t1 in INT_KINDS_TYPES
398+
#:for rank in REDRANKS
399+
#:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp')
400+
module function ${RName}$(x, order, dim, center, mask) result(res)
401+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
402+
integer, intent(in) :: order
403+
integer, intent(in) :: dim
404+
real(dp), intent(in) :: center
405+
logical, intent(in) :: mask${ranksuffix(rank)}$
406+
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
407+
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
427+
428+
end function ${RName}$
429+
#:endfor
430+
#:endfor
431+
432+
276433
#:for k1, t1 in INT_KINDS_TYPES
277434
#:for rank in RANKS
278435
#:set RName = rname("moment_mask",rank, t1, k1, 'dp')

0 commit comments

Comments
 (0)