Skip to content

Addition of variance function in stdlib_experimental_stats #144

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Feb 18, 2020
Merged

Addition of variance function in stdlib_experimental_stats #144

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e966e7b
addition of variance
jvdp1 Feb 5, 2020
044abc5
varaince_dev: update var modules
jvdp1 Feb 5, 2020
d77b6e9
variance_dev: update spec var
jvdp1 Feb 5, 2020
baabfc8
variance_dev: changed ieee_value() as proposed
jvdp1 Feb 7, 2020
9b19154
variance_dev: remove support of complex because it was wrong
jvdp1 Feb 7, 2020
da90a89
variance_dev: addition of variance of complex as (var(real(x)) + var(…
jvdp1 Feb 7, 2020
2a0182a
variance_dev: use fypp to avoid abs in real functions
jvdp1 Feb 7, 2020
01e897c
variance_dev:suggestions by @aradi
jvdp1 Feb 10, 2020
File filter

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
varaince_dev: update var modules
  • Loading branch information
@jvdp1
jvdp1 committed Feb 5, 2020
commit 044abc55653013202c98025fd1ef34d8a23d2adc
14 changes: 14 additions & 0 deletions src/common.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,20 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$
$:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix)
#:enddef

#! Generates Fortran expressions.
#!
#! Args:
#! varname (str): Name of the variable to be used as origin
#! varname1 (str): Name of the variable to be used instead of varname
#! origrank (int): Rank of the original variable
#! dim (int): Index of the used expression varname1
#!
#! Returns:
#! Shape expression enclosed in braces, except for the index dim.
#!
#! E.g., (:, :, :, i, :, :)
#!

#:def rankindice(varname, varname1, origrank, dim)
Copy link
Member

@aradi aradi Feb 10, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should probably rename this macro to have a more descriptive name. If it is only used to select subarrays by reducing the dimension, we could have:

#:def select_subarray(origrank, selectors) #:assert origrank > 0 #:set seldict = dict(selectors) #:call join_lines(joinstr=", ", prefix="(", suffix=")") #:for i in range(1, origrank + 1) $:seldict.get(i, ":") #:endfor #:endcall #:enddef

and use it as

#! -> x(:, i, :) x${select_subarray(3, [(2, 'i')])}$

It could also be used, if we need to reduce more than one rank, e.g.

#! -> x(:, :, i, j) x${select_subarray(4, [(3, 'i'), (4, 'j')])}$

Also the description should be clarified a bit.
Copy link
Member Author

@jvdp1 jvdp1 Feb 10, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Implemented as suggested. The proposed macro is more general and better fit to its aim.
Could you have another review, please?
#:assert origrank > 0
#:if origrank > 0
Expand Down
48 changes: 28 additions & 20 deletions src/stdlib_experimental_stats.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,91 +107,99 @@ module stdlib_experimental_stats

interface var

#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res)
#:set RName = rname("var_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
${t1}$ :: res
end function var_${rank}$_all_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_all_${k1}$_dp(x, mask) result(res)
#:set RName = rname("var_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
real(dp) :: res
end function var_${rank}$_all_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res)
#:set RName = rname("var",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
end function var_${rank}$_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor

#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_${k1}$_dp(x, dim, mask) result(res)
#:set RName = rname("var",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function var_${rank}$_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res)
#:set RName = rname("var_mask_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res
end function var_${rank}$_mask_all_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res)
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res
end function var_${rank}$_mask_all_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res)
#:set RName = rname("var_mask",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
end function var_${rank}$_mask_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res)
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res${reduced_shape('x', rank, 'dim')}$
end function var_${rank}$_mask_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor

Expand Down
63 changes: 37 additions & 26 deletions src/stdlib_experimental_stats_var.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
#:include "common.fypp"

#:set RANKS = range(1, MAXRANK + 1)


#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
submodule (stdlib_experimental_stats) stdlib_experimental_stats_var

use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
Expand All @@ -12,17 +10,18 @@ submodule (stdlib_experimental_stats) stdlib_experimental_stats_var

contains

#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_all_${k1}$_${k1}$(x, mask) result(res)
#:set RName = rname("var_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
${t1}$ :: res

${t1}$ :: n, mean

if (.not.optval(mask, .true.)) then
Copy link
Contributor

@nshaffer nshaffer Feb 11, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a weird idiom to me. Here, I'd prefer the more obvious

if (present(mask)) then if (mask .eqv. .false.) then

But this is a matter of style rather than substance.
Copy link
Member Author

@jvdp1 jvdp1 Feb 13, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hopefully none of both options will be needed in a future standard.
res = ieee_value(res, ieee_quiet_nan)
res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan)
return
end if

Expand All @@ -31,14 +30,15 @@ contains

res = sum((x - mean)**2) / (n - 1._${k1}$)

end function var_${rank}$_all_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_all_${k1}$_dp(x, mask) result(res)
#:set RName = rname("var_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in), optional :: mask
real(dp) :: res
Expand All @@ -55,14 +55,15 @@ contains

res = sum((real(x, dp) - mean)**2) / (n - 1._dp)

end function var_${rank}$_all_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_${k1}$_${k1}$(x, dim, mask) result(res)
#:set RName = rname("var",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
Expand All @@ -73,7 +74,7 @@ contains
${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$

if (.not.optval(mask, .true.)) then
res = ieee_value(res, ieee_quiet_nan)
res = ieee_value(real(res, kind=${k1}$), ieee_quiet_nan)
return
end if

Expand All @@ -92,14 +93,15 @@ contains
end select
res = res / (n - 1._${k1}$)

end function var_${rank}$_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_${k1}$_dp(x, dim, mask) result(res)
#:set RName = rname("var",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in), optional :: mask
Expand Down Expand Up @@ -129,14 +131,15 @@ contains
end select
res = res / (n - 1._dp)

end function var_${rank}$_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_all_${k1}$_${k1}$(x, mask) result(res)
#:set RName = rname("var_mask_all",rank, t1, k1)
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
${t1}$ :: res
Expand All @@ -148,14 +151,15 @@ contains

res = sum((x - mean)**2, mask) / (n - 1._${k1}$)

end function var_${rank}$_mask_all_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_all_${k1}$_dp(x, mask) result(res)
#:set RName = rname("var_mask_all",rank, t1, k1, 'dp')
module function ${RName}$(x, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
logical, intent(in) :: mask${ranksuffix(rank)}$
real(dp) :: res
Expand All @@ -167,14 +171,15 @@ contains

res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp)

end function var_${rank}$_mask_all_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in REAL_KINDS_TYPES
#:for k1, t1 in RC_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_${k1}$_${k1}$(x, dim, mask) result(res)
#:set RName = rname("var_mask",rank, t1, k1)
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
Expand All @@ -192,22 +197,28 @@ contains
mean = sum(x, dim, mask) / n
do i = 1, size(x, dim)
res = res + merge( (x${rankindice(':', 'i', rank, fi )}$ - mean)**2,&
0._${k1}$, mask${rankindice(':', 'i', rank, fi)}$)
#:if t1[0] == 'r'
0._${k1}$,&
#:else
cmplx(0._${k1}$, 0._${k1}$, ${k1}$),&
#:endif
mask${rankindice(':', 'i', rank, fi)}$)
end do
#:endfor
case default
call error_stop("ERROR (mean): wrong dimension")
end select
res = res / (n - 1._${k1}$)

end function var_${rank}$_mask_${k1}$_${k1}$
end function ${RName}$
#:endfor
#:endfor


#:for k1, t1 in INT_KINDS_TYPES
#:for rank in RANKS
module function var_${rank}$_mask_${k1}$_dp(x, dim, mask) result(res)
#:set RName = rname("var_mask",rank, t1, k1, 'dp')
module function ${RName}$(x, dim, mask) result(res)
${t1}$, intent(in) :: x${ranksuffix(rank)}$
integer, intent(in) :: dim
logical, intent(in) :: mask${ranksuffix(rank)}$
Expand All @@ -233,7 +244,7 @@ contains
end select
res = res / (n - 1._dp)

end function var_${rank}$_mask_${k1}$_dp
end function ${RName}$
#:endfor
#:endfor

Expand Down