@@ -11,93 +11,6 @@ submodule (stdlib_stats) stdlib_stats_moment
1111
1212contains
1313
14- #:for k1, t1 in RC_KINDS_TYPES
15- #:for rank in RANKS
16- #:set RName = rname("moment_all",rank, t1, k1)
17- module function ${RName}$(x, order, center, mask) result(res)
18- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
19- integer, intent(in) :: order
20- ${t1}$, intent(in), optional :: center
21- logical, intent(in), optional :: mask
22- ${t1}$ :: res
23-
24- real(${k1}$) :: n
25-
26- if (.not.optval(mask, .true.)) then
27- res = ieee_value(1._${k1}$, ieee_quiet_nan)
28- return
29- end if
30-
31- n = real(size(x, kind = int64), ${k1}$)
32-
33- if (present(center)) then
34- res = sum((x - center)**order) / n
35- else
36- res = sum((x - mean(x))**order) / n
37- end if
38-
39- end function ${RName}$
40- #:endfor
41- #:endfor
42-
43-
44- #:for k1, t1 in INT_KINDS_TYPES
45- #:for rank in RANKS
46- #:set RName = rname("moment_all",rank, t1, k1, 'dp')
47- module function ${RName}$(x, order, center, mask) result(res)
48- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
49- integer, intent(in) :: order
50- real(dp), intent(in), optional :: center
51- logical, intent(in), optional :: mask
52- real(dp) :: res
53-
54- real(dp) :: n
55-
56- if (.not.optval(mask, .true.)) then
57- res = ieee_value(1._dp, ieee_quiet_nan)
58- return
59- end if
60-
61- n = real(size(x, kind = int64), dp)
62-
63- if (present(center)) then
64- res = sum((real(x, dp) - center)**order) / n
65- else
66- res = sum((real(x, dp) - mean(x))**order) / n
67- end if
68-
69- end function ${RName}$
70- #:endfor
71- #:endfor
72-
73-
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- if (.not.optval(mask, .true.)) then
86- res = ieee_value(1._${k1}$, ieee_quiet_nan)
87- return
88- end if
89-
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
95-
96- end function ${RName}$
97- #:endfor
98- #:endfor
99-
100-
10114 #:for k1, t1 in RC_KINDS_TYPES
10215 #:for rank in RANKS
10316 #:set RName = rname("moment",rank, t1, k1)
@@ -146,33 +59,6 @@ contains
14659 #:endfor
14760
14861
149- #:for k1, t1 in INT_KINDS_TYPES
150- #:for rank in REDRANKS
151- #:set RName = rname("moment_scalar",rank, t1, k1, 'dp')
152- module function ${RName}$(x, order, dim, center, mask) result(res)
153- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
154- integer, intent(in) :: order
155- integer, intent(in) :: dim
156- real(dp),intent(in) :: center
157- logical, intent(in), optional :: mask
158- real(dp) :: res${reduced_shape('x', rank, 'dim')}$
159-
160- if (.not.optval(mask, .true.)) then
161- res = ieee_value(1._dp, ieee_quiet_nan)
162- return
163- end if
164-
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
170-
171- end function ${RName}$
172- #:endfor
173- #:endfor
174-
175-
17662 #:for k1, t1 in INT_KINDS_TYPES
17763 #:for rank in RANKS
17864 #:set RName = rname("moment",rank, t1, k1, 'dp')
@@ -221,201 +107,4 @@ contains
221107 #:endfor
222108 #:endfor
223109
224-
225- #:for k1, t1 in RC_KINDS_TYPES
226- #:for rank in RANKS
227- #:set RName = rname("moment_mask_all",rank, t1, k1)
228- module function ${RName}$(x, order, center, mask) result(res)
229- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
230- integer, intent(in) :: order
231- ${t1}$, intent(in), optional :: center
232- logical, intent(in) :: mask${ranksuffix(rank)}$
233- ${t1}$ :: res
234-
235- real(${k1}$) :: n
236-
237- n = real(count(mask, kind = int64), ${k1}$)
238-
239- if (present(center)) then
240- res = sum((x - center)**order, mask) / n
241- else
242- res = sum((x - mean(x, mask))**order, mask) / n
243- end if
244-
245- end function ${RName}$
246- #:endfor
247- #:endfor
248-
249-
250- #:for k1, t1 in INT_KINDS_TYPES
251- #:for rank in RANKS
252- #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp')
253- module function ${RName}$(x, order, center, mask) result(res)
254- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
255- integer, intent(in) :: order
256- real(dp),intent(in), optional :: center
257- logical, intent(in) :: mask${ranksuffix(rank)}$
258- real(dp) :: res
259-
260- real(dp) :: n
261-
262- n = real(count(mask, kind = int64), dp)
263-
264- if (present(center)) then
265- res = sum((real(x, dp) - center)**order, mask) / n
266- else
267- res = sum((real(x, dp) - mean(x,mask))**order, mask) / n
268- end if
269-
270- end function ${RName}$
271- #:endfor
272- #:endfor
273-
274-
275- #:for k1, t1 in RC_KINDS_TYPES
276- #:for rank in REDRANKS
277- #:set RName = rname("moment_mask_scalar",rank, t1, k1)
278- module function ${RName}$(x, order, dim, center, mask) result(res)
279- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
280- integer, intent(in) :: order
281- integer, intent(in) :: dim
282- ${t1}$, intent(in) :: center
283- logical, intent(in) :: mask${ranksuffix(rank)}$
284- ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
285-
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
291-
292- end function ${RName}$
293- #:endfor
294- #:endfor
295-
296-
297- #:for k1, t1 in RC_KINDS_TYPES
298- #:for rank in RANKS
299- #:set RName = rname("moment_mask",rank, t1, k1)
300- module function ${RName}$(x, order, dim, center, mask) result(res)
301- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
302- integer, intent(in) :: order
303- integer, intent(in) :: dim
304- ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
305- logical, intent(in) :: mask${ranksuffix(rank)}$
306- ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
307-
308- integer :: i
309- real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
310- ${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$
311-
312- n = real(count(mask, dim), ${k1}$)
313-
314- res = 0
315- select case(dim)
316- #:for fi in range(1, rank+1)
317- case(${fi}$)
318- if (present(center)) then
319- do i = 1, size(x, ${fi}$)
320- res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ -&
321- center)**order,&
322- #:if t1[0] == 'r'
323- 0._${k1}$,&
324- #:else
325- cmplx(0,0,kind=${k1}$),&
326- #:endif
327- mask${select_subarray(rank, [(fi, 'i')])}$)
328- end do
329- else
330- allocate(mean_, source = mean(x, ${fi}$, mask))
331- do i = 1, size(x, ${fi}$)
332- res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order,&
333- #:if t1[0] == 'r'
334- 0._${k1}$,&
335- #:else
336- cmplx(0,0,kind=${k1}$),&
337- #:endif
338- mask${select_subarray(rank, [(fi, 'i')])}$)
339- end do
340- deallocate(mean_)
341- end if
342- #:endfor
343- case default
344- call error_stop("ERROR (moment): wrong dimension")
345- end select
346- res = res / n
347-
348- end function ${RName}$
349- #:endfor
350- #:endfor
351-
352-
353- #:for k1, t1 in INT_KINDS_TYPES
354- #:for rank in REDRANKS
355- #:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp')
356- module function ${RName}$(x, order, dim, center, mask) result(res)
357- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
358- integer, intent(in) :: order
359- integer, intent(in) :: dim
360- real(dp), intent(in) :: center
361- logical, intent(in) :: mask${ranksuffix(rank)}$
362- real(dp) :: res${reduced_shape('x', rank, 'dim')}$
363-
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
369-
370- end function ${RName}$
371- #:endfor
372- #:endfor
373-
374-
375- #:for k1, t1 in INT_KINDS_TYPES
376- #:for rank in RANKS
377- #:set RName = rname("moment_mask",rank, t1, k1, 'dp')
378- module function ${RName}$(x, order, dim, center, mask) result(res)
379- ${t1}$, intent(in) :: x${ranksuffix(rank)}$
380- integer, intent(in) :: order
381- integer, intent(in) :: dim
382- real(dp), intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
383- logical, intent(in) :: mask${ranksuffix(rank)}$
384- real(dp) :: res${reduced_shape('x', rank, 'dim')}$
385-
386- integer :: i
387- real(dp) :: n${reduced_shape('x', rank, 'dim')}$
388- real(dp), allocatable :: mean_${ranksuffix(rank-1)}$
389-
390- n = real(count(mask, dim), dp)
391-
392- res = 0
393- select case(dim)
394- #:for fi in range(1, rank+1)
395- case(${fi}$)
396- if (present(center)) then
397- do i = 1, size(x, ${fi}$)
398- res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -&
399- center)**order,&
400- 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
401- end do
402- else
403- allocate(mean_, source = mean(x, ${fi}$, mask))
404- do i = 1, size(x, ${fi}$)
405- res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)&
406- **order,&
407- 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
408- end do
409- deallocate(mean_)
410- end if
411- #:endfor
412- case default
413- call error_stop("ERROR (moment): wrong dimension")
414- end select
415- res = res / n
416-
417- end function ${RName}$
418- #:endfor
419- #:endfor
420-
421110end submodule
0 commit comments