Skip to content

Commit 19c1f8e

Browse files
jvdp1jvdp1
committed
corr_dev: implemented some correlation functions
1 parent ec16e9c commit 19c1f8e

File tree

1 file changed

+137
-131
lines changed

1 file changed

+137
-131
lines changed

‎src/stdlib_experimental_stats_corr.fypp

Lines changed: 137 additions & 131 deletions
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,10 @@ contains
1111

1212
#:for k1, t1 in RC_KINDS_TYPES
1313
#:set RName = rname("corr",1, t1, k1)
14-
module function ${RName}$(x, dim, mask, corrected) result(res)
14+
module function ${RName}$(x, dim, mask) result(res)
1515
${t1}$, intent(in) :: x(:)
1616
integer, intent(in) :: dim
1717
logical, intent(in), optional :: mask
18-
logical, intent(in), optional :: corrected
1918
real(${k1}$) :: res
2019

2120
if (.not.optval(mask, .true.)) then
@@ -31,11 +30,10 @@ contains
3130

3231
#:for k1, t1 in INT_KINDS_TYPES
3332
#:set RName = rname("corr",1, t1, k1, 'dp')
34-
module function ${RName}$(x, dim, mask, corrected) result(res)
33+
module function ${RName}$(x, dim, mask) result(res)
3534
${t1}$, intent(in) :: x(:)
3635
integer, intent(in) :: dim
3736
logical, intent(in), optional :: mask
38-
logical, intent(in), optional :: corrected
3937
real(dp) :: res
4038

4139
if (.not.optval(mask, .true.)) then
@@ -51,11 +49,10 @@ contains
5149

5250
#:for k1, t1 in RC_KINDS_TYPES
5351
#:set RName = rname("corr_mask",1, t1, k1)
54-
module function ${RName}$(x, dim, mask, corrected) result(res)
52+
module function ${RName}$(x, dim, mask) result(res)
5553
${t1}$, intent(in) :: x(:)
5654
integer, intent(in) :: dim
5755
logical, intent(in) :: mask(:)
58-
logical, intent(in), optional :: corrected
5956
real(${k1}$) :: res
6057

6158
res = 1
@@ -66,11 +63,10 @@ contains
6663

6764
#:for k1, t1 in INT_KINDS_TYPES
6865
#:set RName = rname("corr_mask",1, t1, k1, 'dp')
69-
module function ${RName}$(x, dim, mask, corrected) result(res)
66+
module function ${RName}$(x, dim, mask) result(res)
7067
${t1}$, intent(in) :: x(:)
7168
integer, intent(in) :: dim
7269
logical, intent(in) :: mask(:)
73-
logical, intent(in), optional :: corrected
7470
real(dp) :: res
7571

7672
res = 1
@@ -81,11 +77,10 @@ contains
8177

8278
#:for k1, t1 in RC_KINDS_TYPES
8379
#:set RName = rname("corr",2, t1, k1)
84-
module function ${RName}$(x, dim, mask, corrected) result(res)
80+
module function ${RName}$(x, dim, mask) result(res)
8581
${t1}$, intent(in) :: x(:, :)
8682
integer, intent(in) :: dim
8783
logical, intent(in), optional :: mask
88-
logical, intent(in), optional :: corrected
8984
${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
9085
, merge(size(x, 1), size(x, 2), mask = 1<dim))
9186

@@ -121,19 +116,24 @@ contains
121116
case default
122117
call error_stop("ERROR (corr): wrong dimension")
123118
end select
124-
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
119+
120+
mean_ = 1 / sqrt(diag(res))
121+
do i = 1, size(res, 1)
122+
do j = 1, size(res, 2)
123+
res(j, i) = res(j, i) * mean_(i) * mean_(j)
124+
end do
125+
end do
125126

126127
end function ${RName}$
127128
#:endfor
128129

129130

130131
#:for k1, t1 in INT_KINDS_TYPES
131132
#:set RName = rname("corr",2, t1, k1, 'dp')
132-
module function ${RName}$(x, dim, mask, corrected) result(res)
133+
module function ${RName}$(x, dim, mask) result(res)
133134
${t1}$, intent(in) :: x(:, :)
134135
integer, intent(in) :: dim
135136
logical, intent(in), optional :: mask
136-
logical, intent(in), optional :: corrected
137137
real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
138138
, merge(size(x, 1), size(x, 2), mask = 1<dim))
139139

@@ -161,130 +161,136 @@ contains
161161
case default
162162
call error_stop("ERROR (corr): wrong dimension")
163163
end select
164-
res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
165164

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
165+
mean_ = 1 / sqrt(diag(res))
166+
do i = 1, size(res, 1)
167+
do j = 1, size(res, 2)
168+
res(j, i) = res(j, i) * mean_(i) * mean_(j)
169+
end do
170+
end do
233171

234172
end function ${RName}$
235173
#:endfor
236174

237175

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
176+
! #:for k1, t1 in RC_KINDS_TYPES
177+
! #:set RName = rname("corr_mask",2, t1, k1)
178+
! module function ${RName}$(x, dim, mask, corrected) result(res)
179+
! ${t1}$, intent(in) :: x(:, :)
180+
! integer, intent(in) :: dim
181+
! logical, intent(in) :: mask(:,:)
182+
! logical, intent(in), optional :: corrected
183+
! ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
184+
! , merge(size(x, 1), size(x, 2), mask = 1<dim))
185+
!
186+
! integer :: i, j, n
187+
! ${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
188+
! ${t1}$ :: center(size(x, 1),size(x, 2))
189+
!
190+
! mean_ = mean(x, dim, mask = mask)
191+
! select case(dim)
192+
! case(1)
193+
! do i = 1, size(x, 1)
194+
! center(i, :) = merge( x(i, :) - mean_,&
195+
! #:if t1[0] == 'r'
196+
! 0._${k1}$,&
197+
! #:else
198+
! cmplx(0,0,kind=${k1}$),&
199+
! #:endif
200+
! mask(i, :))
201+
! end do
202+
! #:if t1[0] == 'r'
203+
! res = matmul( transpose(center), center)
204+
! #:else
205+
! res = matmul( transpose(conjg(center)), center)
206+
! #:endif
207+
! do j = 1, size(res, 2)
208+
! do i = 1, size(res, 1)
209+
! n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
210+
! res(i, j) = res(i, j) / (n - merge(1, 0,&
211+
! optval(corrected, .true.) .and. n > 0))
212+
! end do
213+
! end do
214+
! case(2)
215+
! do i = 1, size(x, 2)
216+
! center(:, i) = merge( x(:, i) - mean_,&
217+
! #:if t1[0] == 'r'
218+
! 0._${k1}$,&
219+
! #:else
220+
! cmplx(0,0,kind=${k1}$),&
221+
! #:endif
222+
! mask(:, i))
223+
! end do
224+
! #:if t1[0] == 'r'
225+
! res = matmul( center, transpose(center))
226+
! #:else
227+
! res = matmul( center, transpose(conjg(center)))
228+
! #:endif
229+
! do j = 1, size(res, 2)
230+
! do i = 1, size(res, 1)
231+
! n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
232+
! res(i, j) = res(i, j) / (n - merge(1, 0,&
233+
! optval(corrected, .true.) .and. n > 0))
234+
! end do
235+
! end do
236+
! case default
237+
! call error_stop("ERROR (corr): wrong dimension")
238+
! end select
239+
!
240+
! end function ${RName}$
241+
! #:endfor
242+
!
243+
!
244+
! #:for k1, t1 in INT_KINDS_TYPES
245+
! #:set RName = rname("corr_mask",2, t1, k1, 'dp')
246+
! module function ${RName}$(x, dim, mask, corrected) result(res)
247+
! ${t1}$, intent(in) :: x(:, :)
248+
! integer, intent(in) :: dim
249+
! logical, intent(in) :: mask(:,:)
250+
! logical, intent(in), optional :: corrected
251+
! real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
252+
! , merge(size(x, 1), size(x, 2), mask = 1<dim))
253+
!
254+
! integer :: i, j, n
255+
! real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
256+
! real(dp) :: center(size(x, 1),size(x, 2))
257+
!
258+
! mean_ = mean(x, dim, mask = mask)
259+
! select case(dim)
260+
! case(1)
261+
! do i = 1, size(x, 1)
262+
! center(i, :) = merge( x(i, :) - mean_,&
263+
! 0._dp,&
264+
! mask(i, :))
265+
! end do
266+
! res = matmul( transpose(center), center)
267+
! do j = 1, size(res, 2)
268+
! do i = 1, size(res, 1)
269+
! n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
270+
! res(i, j) = res(i, j) / (n - merge(1, 0,&
271+
! optval(corrected, .true.) .and. n > 0))
272+
! end do
273+
! end do
274+
! case(2)
275+
! do i = 1, size(x, 2)
276+
! center(:, i) = merge( x(:, i) - mean_,&
277+
! 0._dp,&
278+
! mask(:, i))
279+
! end do
280+
! res = matmul( center, transpose(center))
281+
! do j = 1, size(res, 2)
282+
! do i = 1, size(res, 1)
283+
! n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
284+
! res(i, j) = res(i, j) / (n - merge(1, 0,&
285+
! optval(corrected, .true.) .and. n > 0))
286+
! end do
287+
! end do
288+
! case default
289+
! call error_stop("ERROR (corr): wrong dimension")
290+
! end select
291+
!
292+
! end function ${RName}$
293+
! #:endfor
288294

289295

290296
end submodule

0 commit comments

Comments
 (0)