Skip to content

Commit ec16e9c

Browse files
jvdp1jvdp1
committed
corr_dev: init
1 parent 9e377fe commit ec16e9c

File tree

1 file changed

+290
-0
lines changed

1 file changed

+290
-0
lines changed

‎src/stdlib_experimental_stats_corr.fypp

Lines changed: 290 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,290 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
submodule (stdlib_experimental_stats) stdlib_experimental_stats_corr
4+
5+
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
6+
use stdlib_experimental_error, only: error_stop
7+
use stdlib_experimental_optval, only: optval
8+
implicit none
9+
10+
contains
11+
12+
#:for k1, t1 in RC_KINDS_TYPES
13+
#:set RName = rname("corr",1, t1, k1)
14+
module function ${RName}$(x, dim, mask, corrected) result(res)
15+
${t1}$, intent(in) :: x(:)
16+
integer, intent(in) :: dim
17+
logical, intent(in), optional :: mask
18+
logical, intent(in), optional :: corrected
19+
real(${k1}$) :: res
20+
21+
if (.not.optval(mask, .true.)) then
22+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
23+
return
24+
end if
25+
26+
res = 1
27+
28+
end function ${RName}$
29+
#:endfor
30+
31+
32+
#:for k1, t1 in INT_KINDS_TYPES
33+
#:set RName = rname("corr",1, t1, k1, 'dp')
34+
module function ${RName}$(x, dim, mask, corrected) result(res)
35+
${t1}$, intent(in) :: x(:)
36+
integer, intent(in) :: dim
37+
logical, intent(in), optional :: mask
38+
logical, intent(in), optional :: corrected
39+
real(dp) :: res
40+
41+
if (.not.optval(mask, .true.)) then
42+
res = ieee_value(1._dp, ieee_quiet_nan)
43+
return
44+
end if
45+
46+
res = 1
47+
48+
end function ${RName}$
49+
#:endfor
50+
51+
52+
#:for k1, t1 in RC_KINDS_TYPES
53+
#:set RName = rname("corr_mask",1, t1, k1)
54+
module function ${RName}$(x, dim, mask, corrected) result(res)
55+
${t1}$, intent(in) :: x(:)
56+
integer, intent(in) :: dim
57+
logical, intent(in) :: mask(:)
58+
logical, intent(in), optional :: corrected
59+
real(${k1}$) :: res
60+
61+
res = 1
62+
63+
end function ${RName}$
64+
#:endfor
65+
66+
67+
#:for k1, t1 in INT_KINDS_TYPES
68+
#:set RName = rname("corr_mask",1, t1, k1, 'dp')
69+
module function ${RName}$(x, dim, mask, corrected) result(res)
70+
${t1}$, intent(in) :: x(:)
71+
integer, intent(in) :: dim
72+
logical, intent(in) :: mask(:)
73+
logical, intent(in), optional :: corrected
74+
real(dp) :: res
75+
76+
res = 1
77+
78+
end function ${RName}$
79+
#:endfor
80+
81+
82+
#:for k1, t1 in RC_KINDS_TYPES
83+
#:set RName = rname("corr",2, t1, k1)
84+
module function ${RName}$(x, dim, mask, corrected) result(res)
85+
${t1}$, intent(in) :: x(:, :)
86+
integer, intent(in) :: dim
87+
logical, intent(in), optional :: mask
88+
logical, intent(in), optional :: corrected
89+
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
90+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
91+
92+
integer :: i
93+
${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
94+
${t1}$ :: center(size(x, 1),size(x, 2))
95+
96+
if (.not.optval(mask, .true.)) then
97+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
98+
return
99+
end if
100+
101+
mean_ = mean(x, dim)
102+
select case(dim)
103+
case(1)
104+
do i = 1, size(x, 1)
105+
center(i, :) = x(i, :) - mean_
106+
end do
107+
#:if t1[0] == 'r'
108+
res = matmul( transpose(center), center)
109+
#:else
110+
res = matmul( transpose(conjg(center)), center)
111+
#:endif
112+
case(2)
113+
do i = 1, size(x, 2)
114+
center(:, i) = x(:, i) - mean_
115+
end do
116+
#:if t1[0] == 'r'
117+
res = matmul( center, transpose(center))
118+
#:else
119+
res = matmul( center, transpose(conjg(center)))
120+
#:endif
121+
case default
122+
call error_stop("ERROR (corr): wrong dimension")
123+
end select
124+
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
125+
126+
end function ${RName}$
127+
#:endfor
128+
129+
130+
#:for k1, t1 in INT_KINDS_TYPES
131+
#:set RName = rname("corr",2, t1, k1, 'dp')
132+
module function ${RName}$(x, dim, mask, corrected) result(res)
133+
${t1}$, intent(in) :: x(:, :)
134+
integer, intent(in) :: dim
135+
logical, intent(in), optional :: mask
136+
logical, intent(in), optional :: corrected
137+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
138+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
139+
140+
integer :: i
141+
real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
142+
real(dp) :: center(size(x, 1),size(x, 2))
143+
144+
if (.not.optval(mask, .true.)) then
145+
res = ieee_value(1._dp, ieee_quiet_nan)
146+
return
147+
end if
148+
149+
mean_ = mean(x, dim)
150+
select case(dim)
151+
case(1)
152+
do i = 1, size(x, 1)
153+
center(i, :) = real(x(i, :), dp) - mean_
154+
end do
155+
res = matmul( transpose(center), center)
156+
case(2)
157+
do i = 1, size(x, 2)
158+
center(:, i) = real(x(:, i), dp) - mean_
159+
end do
160+
res = matmul( center, transpose(center))
161+
case default
162+
call error_stop("ERROR (corr): wrong dimension")
163+
end select
164+
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
165+
166+
end function ${RName}$
167+
#:endfor
168+
169+
170+
#:for k1, t1 in RC_KINDS_TYPES
171+
#:set RName = rname("corr_mask",2, t1, k1)
172+
module function ${RName}$(x, dim, mask, corrected) result(res)
173+
${t1}$, intent(in) :: x(:, :)
174+
integer, intent(in) :: dim
175+
logical, intent(in) :: mask(:,:)
176+
logical, intent(in), optional :: corrected
177+
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
178+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
179+
180+
integer :: i, j, n
181+
${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
182+
${t1}$ :: center(size(x, 1),size(x, 2))
183+
184+
mean_ = mean(x, dim, mask = mask)
185+
select case(dim)
186+
case(1)
187+
do i = 1, size(x, 1)
188+
center(i, :) = merge( x(i, :) - mean_,&
189+
#:if t1[0] == 'r'
190+
0._${k1}$,&
191+
#:else
192+
cmplx(0,0,kind=${k1}$),&
193+
#:endif
194+
mask(i, :))
195+
end do
196+
#:if t1[0] == 'r'
197+
res = matmul( transpose(center), center)
198+
#:else
199+
res = matmul( transpose(conjg(center)), center)
200+
#:endif
201+
do j = 1, size(res, 2)
202+
do i = 1, size(res, 1)
203+
n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
204+
res(i, j) = res(i, j) / (n - merge(1, 0,&
205+
optval(corrected, .true.) .and. n > 0))
206+
end do
207+
end do
208+
case(2)
209+
do i = 1, size(x, 2)
210+
center(:, i) = merge( x(:, i) - mean_,&
211+
#:if t1[0] == 'r'
212+
0._${k1}$,&
213+
#:else
214+
cmplx(0,0,kind=${k1}$),&
215+
#:endif
216+
mask(:, i))
217+
end do
218+
#:if t1[0] == 'r'
219+
res = matmul( center, transpose(center))
220+
#:else
221+
res = matmul( center, transpose(conjg(center)))
222+
#:endif
223+
do j = 1, size(res, 2)
224+
do i = 1, size(res, 1)
225+
n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
226+
res(i, j) = res(i, j) / (n - merge(1, 0,&
227+
optval(corrected, .true.) .and. n > 0))
228+
end do
229+
end do
230+
case default
231+
call error_stop("ERROR (corr): wrong dimension")
232+
end select
233+
234+
end function ${RName}$
235+
#:endfor
236+
237+
238+
#:for k1, t1 in INT_KINDS_TYPES
239+
#:set RName = rname("corr_mask",2, t1, k1, 'dp')
240+
module function ${RName}$(x, dim, mask, corrected) result(res)
241+
${t1}$, intent(in) :: x(:, :)
242+
integer, intent(in) :: dim
243+
logical, intent(in) :: mask(:,:)
244+
logical, intent(in), optional :: corrected
245+
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
246+
, merge(size(x, 1), size(x, 2), mask = 1<dim))
247+
248+
integer :: i, j, n
249+
real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
250+
real(dp) :: center(size(x, 1),size(x, 2))
251+
252+
mean_ = mean(x, dim, mask = mask)
253+
select case(dim)
254+
case(1)
255+
do i = 1, size(x, 1)
256+
center(i, :) = merge( x(i, :) - mean_,&
257+
0._dp,&
258+
mask(i, :))
259+
end do
260+
res = matmul( transpose(center), center)
261+
do j = 1, size(res, 2)
262+
do i = 1, size(res, 1)
263+
n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
264+
res(i, j) = res(i, j) / (n - merge(1, 0,&
265+
optval(corrected, .true.) .and. n > 0))
266+
end do
267+
end do
268+
case(2)
269+
do i = 1, size(x, 2)
270+
center(:, i) = merge( x(:, i) - mean_,&
271+
0._dp,&
272+
mask(:, i))
273+
end do
274+
res = matmul( center, transpose(center))
275+
do j = 1, size(res, 2)
276+
do i = 1, size(res, 1)
277+
n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
278+
res(i, j) = res(i, j) / (n - merge(1, 0,&
279+
optval(corrected, .true.) .and. n > 0))
280+
end do
281+
end do
282+
case default
283+
call error_stop("ERROR (corr): wrong dimension")
284+
end select
285+
286+
end function ${RName}$
287+
#:endfor
288+
289+
290+
end submodule

0 commit comments

Comments
 (0)