Skip to content

Commit da90a89

Browse files
jvdp1jvdp1
committed
variance_dev: addition of variance of complex as (var(real(x)) + var(aimag(x))
1 parent 9b19154 commit da90a89

File tree

3 files changed

+278
-7
lines changed

3 files changed

+278
-7
lines changed

‎src/stdlib_experimental_stats.fypp‎

Lines changed: 48 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ module stdlib_experimental_stats
107107

108108
interface var
109109

110-
#:for k1, t1 in RC_KINDS_TYPES
110+
#:for k1, t1 in REAL_KINDS_TYPES
111111
#:for rank in RANKS
112112
#:set RName = rname("var_all",rank, t1, k1)
113113
module function ${RName}$(x, mask) result(res)
@@ -118,6 +118,17 @@ module stdlib_experimental_stats
118118
#:endfor
119119
#:endfor
120120

121+
#:for k1, t1 in CMPLX_KINDS_TYPES
122+
#:for rank in RANKS
123+
#:set RName = rname("var_all",rank, t1, k1, 'r'+k1)
124+
module function ${RName}$(x, mask) result(res)
125+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
126+
logical, intent(in), optional :: mask
127+
real(${k1}$) :: res
128+
end function ${RName}$
129+
#:endfor
130+
#:endfor
131+
121132
#:for k1, t1 in INT_KINDS_TYPES
122133
#:for rank in RANKS
123134
#:set RName = rname("var_all",rank, t1, k1, 'dp')
@@ -129,7 +140,7 @@ module stdlib_experimental_stats
129140
#:endfor
130141
#:endfor
131142

132-
#:for k1, t1 in RC_KINDS_TYPES
143+
#:for k1, t1 in REAL_KINDS_TYPES
133144
#:for rank in RANKS
134145
#:set RName = rname("var",rank, t1, k1)
135146
module function ${RName}$(x, dim, mask) result(res)
@@ -141,6 +152,18 @@ module stdlib_experimental_stats
141152
#:endfor
142153
#:endfor
143154

155+
#:for k1, t1 in CMPLX_KINDS_TYPES
156+
#:for rank in RANKS
157+
#:set RName = rname("var",rank, t1, k1, 'r'+k1)
158+
module function ${RName}$(x, dim, mask) result(res)
159+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
160+
integer, intent(in) :: dim
161+
logical, intent(in), optional :: mask
162+
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
163+
end function ${RName}$
164+
#:endfor
165+
#:endfor
166+
144167
#:for k1, t1 in INT_KINDS_TYPES
145168
#:for rank in RANKS
146169
#:set RName = rname("var",rank, t1, k1, 'dp')
@@ -153,8 +176,7 @@ module stdlib_experimental_stats
153176
#:endfor
154177
#:endfor
155178

156-
157-
#:for k1, t1 in RC_KINDS_TYPES
179+
#:for k1, t1 in REAL_KINDS_TYPES
158180
#:for rank in RANKS
159181
#:set RName = rname("var_mask_all",rank, t1, k1)
160182
module function ${RName}$(x, mask) result(res)
@@ -165,6 +187,16 @@ module stdlib_experimental_stats
165187
#:endfor
166188
#:endfor
167189

190+
#:for k1, t1 in CMPLX_KINDS_TYPES
191+
#:for rank in RANKS
192+
#:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1)
193+
module function ${RName}$(x, mask) result(res)
194+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
195+
logical, intent(in) :: mask${ranksuffix(rank)}$
196+
real(${k1}$) :: res
197+
end function ${RName}$
198+
#:endfor
199+
#:endfor
168200

169201
#:for k1, t1 in INT_KINDS_TYPES
170202
#:for rank in RANKS
@@ -177,8 +209,7 @@ module stdlib_experimental_stats
177209
#:endfor
178210
#:endfor
179211

180-
181-
#:for k1, t1 in RC_KINDS_TYPES
212+
#:for k1, t1 in REAL_KINDS_TYPES
182213
#:for rank in RANKS
183214
#:set RName = rname("var_mask",rank, t1, k1)
184215
module function ${RName}$(x, dim, mask) result(res)
@@ -190,6 +221,17 @@ module stdlib_experimental_stats
190221
#:endfor
191222
#:endfor
192223

224+
#:for k1, t1 in CMPLX_KINDS_TYPES
225+
#:for rank in RANKS
226+
#:set RName = rname("var_mask",rank, t1, k1, 'r'+k1)
227+
module function ${RName}$(x, dim, mask) result(res)
228+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
229+
integer, intent(in) :: dim
230+
logical, intent(in) :: mask${ranksuffix(rank)}$
231+
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
232+
end function ${RName}$
233+
#:endfor
234+
#:endfor
193235

194236
#:for k1, t1 in INT_KINDS_TYPES
195237
#:for rank in RANKS

‎src/stdlib_experimental_stats_var.fypp‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,36 @@ contains
3434
#:endfor
3535

3636

37+
#:for k1, t1 in CMPLX_KINDS_TYPES
38+
#:for rank in RANKS
39+
#:set RName = rname("var_all",rank, t1, k1, 'r'+k1)
40+
module function ${RName}$(x, mask) result(res)
41+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
42+
logical, intent(in), optional :: mask
43+
real(${k1}$) :: res
44+
45+
real(${k1}$) :: n, mean
46+
47+
if (.not.optval(mask, .true.)) then
48+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
49+
return
50+
end if
51+
52+
n = real(size(x, kind = int64), ${k1}$)
53+
54+
!real part
55+
mean = sum(real(x)) / n
56+
res = sum((real(x) - mean)**2)
57+
58+
!imaginary part
59+
mean = sum(aimag(x)) / n
60+
res = (res + sum((aimag(x) - mean)**2)) / (n - 1._${k1}$)
61+
62+
end function ${RName}$
63+
#:endfor
64+
#:endfor
65+
66+
3767
#:for k1, t1 in INT_KINDS_TYPES
3868
#:for rank in RANKS
3969
#:set RName = rname("var_all",rank, t1, k1, 'dp')
@@ -97,6 +127,50 @@ contains
97127
#:endfor
98128

99129

130+
#:for k1, t1 in CMPLX_KINDS_TYPES
131+
#:for rank in RANKS
132+
#:set RName = rname("var",rank, t1, k1, 'r'+k1)
133+
module function ${RName}$(x, dim, mask) result(res)
134+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
135+
integer, intent(in) :: dim
136+
logical, intent(in), optional :: mask
137+
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
138+
139+
integer :: i
140+
real(${k1}$) :: n
141+
real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$
142+
143+
if (.not.optval(mask, .true.)) then
144+
res = ieee_value(1._${k1}$, ieee_quiet_nan)
145+
return
146+
end if
147+
148+
res = 0._${k1}$
149+
select case(dim)
150+
#:for fi in range(1, rank+1)
151+
case(${fi}$)
152+
n = real(size(x, dim), ${k1}$)
153+
!real part
154+
mean = sum(real(x), dim) / n
155+
do i = 1, size(x, dim)
156+
res = res + (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2
157+
end do
158+
!imaginary part
159+
mean = sum(aimag(x), dim) / n
160+
do i = 1, size(x, dim)
161+
res = res + (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2
162+
end do
163+
#:endfor
164+
case default
165+
call error_stop("ERROR (mean): wrong dimension")
166+
end select
167+
res = res / (n - 1._${k1}$)
168+
169+
end function ${RName}$
170+
#:endfor
171+
#:endfor
172+
173+
100174
#:for k1, t1 in INT_KINDS_TYPES
101175
#:for rank in RANKS
102176
#:set RName = rname("var",rank, t1, k1, 'dp')
@@ -155,6 +229,31 @@ contains
155229
#:endfor
156230

157231

232+
#:for k1, t1 in CMPLX_KINDS_TYPES
233+
#:for rank in RANKS
234+
#:set RName = rname("var_mask_all",rank, t1, k1, 'r'+k1)
235+
module function ${RName}$(x, mask) result(res)
236+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
237+
logical, intent(in) :: mask${ranksuffix(rank)}$
238+
real(${k1}$) :: res
239+
240+
real(${k1}$) :: n, mean
241+
242+
n = real(count(mask, kind = int64), ${k1}$)
243+
244+
!real part
245+
mean = sum(real(x), mask) / n
246+
res = sum((real(x) - mean)**2, mask)
247+
248+
!imaginary part
249+
mean = sum(aimag(x), mask) / n
250+
res = (res + sum((aimag(x) - mean)**2, mask)) / (n - 1._${k1}$)
251+
252+
end function ${RName}$
253+
#:endfor
254+
#:endfor
255+
256+
158257
#:for k1, t1 in INT_KINDS_TYPES
159258
#:for rank in RANKS
160259
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
@@ -210,6 +309,49 @@ contains
210309
#:endfor
211310

212311

312+
#:for k1, t1 in CMPLX_KINDS_TYPES
313+
#:for rank in RANKS
314+
#:set RName = rname("var_mask",rank, t1, k1, 'r'+k1)
315+
module function ${RName}$(x, dim, mask) result(res)
316+
${t1}$, intent(in) :: x${ranksuffix(rank)}$
317+
integer, intent(in) :: dim
318+
logical, intent(in) :: mask${ranksuffix(rank)}$
319+
real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$
320+
321+
integer :: i
322+
real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
323+
real(${k1}$) :: mean${reduced_shape('x', rank, 'dim')}$
324+
325+
res = 0._${k1}$
326+
select case(dim)
327+
#:for fi in range(1, rank+1)
328+
case(${fi}$)
329+
n = real(count(mask, dim), ${k1}$)
330+
!real part
331+
mean = sum(real(x), dim, mask) / n
332+
do i = 1, size(x, dim)
333+
res = res + merge( (real(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,&
334+
0._${k1}$,&
335+
mask${rankindice(':', 'i', rank, fi)}$)
336+
end do
337+
!imaginary part
338+
mean = sum(aimag(x), dim, mask) / n
339+
do i = 1, size(x, dim)
340+
res = res + merge( (aimag(x${rankindice(':', 'i', rank, fi )}$) - mean)**2,&
341+
0._${k1}$,&
342+
mask${rankindice(':', 'i', rank, fi)}$)
343+
end do
344+
#:endfor
345+
case default
346+
call error_stop("ERROR (mean): wrong dimension")
347+
end select
348+
res = res / (n - 1._${k1}$)
349+
350+
end function ${RName}$
351+
#:endfor
352+
#:endfor
353+
354+
213355
#:for k1, t1 in INT_KINDS_TYPES
214356
#:for rank in RANKS
215357
#:set RName = rname("var_mask",rank, t1, k1, 'dp')

‎src/tests/stats/test_var.f90‎

Lines changed: 88 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,10 +20,25 @@ program test_var
2020

2121
real(sp), allocatable :: s(:, :), s3(:, :, :)
2222
real(dp), allocatable :: d3(:, :, :)
23-
real(dp) :: d(4, 3)= reshape([1._dp, 3._dp, 5._dp, 7._dp,&
23+
real(dp) :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 7._dp,&
2424
2._dp, 4._dp, 6._dp, 8._dp,&
2525
9._dp, 10._dp, 11._dp, 12._dp], [4, 3])
2626

27+
28+
complex(sp) :: cs1(5) = [ cmplx(0.57706_sp, 0.00000_sp),&
29+
cmplx(0.00000_sp, 1.44065_sp),&
30+
cmplx(1.26401_sp, 0.00000_sp),&
31+
cmplx(0.00000_sp, 0.88833_sp),&
32+
cmplx(1.14352_sp, 0.00000_sp)]
33+
complex(dp) :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp),&
34+
cmplx(0.00000_dp, 1.44065_dp),&
35+
cmplx(1.26401_dp, 0.00000_dp),&
36+
cmplx(0.00000_dp, 0.88833_dp),&
37+
cmplx(1.14352_dp, 0.00000_dp)]
38+
complex(sp) :: cs(5,3)
39+
complex(dp) :: cd(5,3)
40+
41+
2742
!sp
2843
!1dim
2944
print*,' test_sp_1dim'
@@ -388,4 +403,76 @@ program test_var
388403
[size(i643,1), size(i643,2)] ))&
389404
< dptol ))
390405

406+
!csp
407+
!1dim
408+
print*,' test_csp_1dim'
409+
call assert( abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol)
410+
call assert( abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol)
411+
412+
print*,' test_csp_1dim_mask'
413+
call assert( isnan(var(cs1, .false.)))
414+
call assert( isnan(var(cs1, 1, .false.)))
415+
416+
print*,' test_sp_1dim_mask_array'
417+
call assert( abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol)
418+
call assert( abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol)
419+
420+
!2dim
421+
cs(:,1) = cs1
422+
cs(:,2) = cs1*3_sp
423+
cs(:,3) = cs1*1.5_sp
424+
425+
print*,' test_sp_2dim'
426+
print*,var(cs),(var(real(cs)) + var(aimag(cs)))
427+
print*,var(cs)-(var(real(cs)) + var(aimag(cs))),sptol
428+
call assert( abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol)
429+
call assert( all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol))
430+
call assert( all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol))
431+
432+
print*,' test_sp_2dim_mask'
433+
call assert( isnan(var(cs, .false.)))
434+
call assert( any(isnan(var(cs, 1, .false.))))
435+
call assert( any(isnan(var(cs, 2, .false.))))
436+
437+
print*,' test_sp_2dim_mask_array'
438+
call assert( abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol)
439+
call assert( all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol))
440+
call assert( all( abs( var(cs, 2, aimag(cs) == 0) - var(real(cs), 2, aimag(cs) == 0)) < sptol))
441+
442+
!cdp
443+
!1dim
444+
print*,' test_cdp_1dim'
445+
call assert( abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol)
446+
call assert( abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol)
447+
448+
print*,' test_cdp_1dim_mask'
449+
call assert( isnan(var(cd1, .false.)))
450+
call assert( isnan(var(cd1, 1, .false.)))
451+
452+
print*,' test_sp_1dim_mask_array'
453+
call assert( abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol)
454+
call assert( abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol)
455+
456+
!2dim
457+
cd(:,1) = cd1
458+
cd(:,2) = cd1*3_sp
459+
cd(:,3) = cd1*1.5_sp
460+
461+
print*,' test_sp_2dim'
462+
print*,var(cd),(var(real(cd)) + var(aimag(cd)))
463+
print*,var(cd)-(var(real(cd)) + var(aimag(cd))),dptol
464+
call assert( abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol)
465+
call assert( all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol))
466+
call assert( all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol))
467+
468+
print*,' test_sp_2dim_mask'
469+
call assert( isnan(var(cd, .false.)))
470+
call assert( any(isnan(var(cd, 1, .false.))))
471+
call assert( any(isnan(var(cd, 2, .false.))))
472+
473+
print*,' test_sp_2dim_mask_array'
474+
call assert( abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol)
475+
call assert( all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol))
476+
call assert( all( abs( var(cd, 2, aimag(cd) == 0) - var(real(cd), 2, aimag(cd) == 0)) < dptol))
477+
391478
end program

0 commit comments

Comments
 (0)