Skip to content

convert pdf_norm and cdf_norm to pure while improving scale check #679

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
Mar 4, 2023
Merged

convert pdf_norm and cdf_norm to pure while improving scale check #679

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d4e6c67
convert proc to pure, comment out scale check
HugoMVale Sep 14, 2022
f60d75a
switch back rvs_norm to impure
HugoMVale Sep 14, 2022
a2b8932
remove comments and prettify
HugoMVale Sep 16, 2022
ecfd440
added protection against scale<=0
HugoMVale Sep 16, 2022
acb8fad
fix nan implementation
HugoMVale Sep 16, 2022
a656a51
add ieee_value
HugoMVale Sep 16, 2022
f63afaf
fprettify -i 4
HugoMVale Sep 21, 2022
a2e8923
Merge remote-tracking branch 'upstream/master' into stats_dist_normal…
milancurcic Mar 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 34 additions & 47 deletions src/stdlib_stats_distribution_normal.fypp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
module stdlib_stats_distribution_normal
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
use stdlib_kinds, only: sp, dp, xdp, qp, int32
use stdlib_error, only: error_stop
use stdlib_random, only: dist_rand
Expand All @@ -18,8 +19,6 @@ module stdlib_stats_distribution_normal
public :: pdf_normal
public :: cdf_normal



interface rvs_normal
!! version: experimental
!!
Expand All @@ -38,8 +37,6 @@ module stdlib_stats_distribution_normal
#:endfor
end interface rvs_normal



interface pdf_normal
!! version: experimental
!!
Expand All @@ -52,8 +49,6 @@ module stdlib_stats_distribution_normal
#:endfor
end interface pdf_normal



interface cdf_normal
!! version: experimental
!!
Expand All @@ -66,13 +61,9 @@ module stdlib_stats_distribution_normal
#:endfor
end interface cdf_normal





contains

subroutine zigset
impure subroutine zigset
Copy link
Member

Choose a reason for hiding this comment

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

This impure attribute is redundant. I would suggest to remove it in this case and other occurrences.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That deserves a bit of explanation. IMO, the attribute is not redundant from a "documentation" perspective.
A function without an explicit attribute will be treated as impure, but it may in reality be pure or impure. The only way to find out it is to look inside the code. Adding the impure statement there is meant to explicitly indicate that that function is really impure (it changes global variables).

Copy link
Member

Choose a reason for hiding this comment

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

I think it's okay for impure to stay for documentation purposes.

! Marsaglia & Tsang generator for random normals & random exponentials.
! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public
! domain (https://jblevins.org/mirror/amiller/)
Expand Down Expand Up @@ -109,12 +100,10 @@ contains
zig_norm_initialized = .true.
end subroutine zigset



#:for k1, t1 in REAL_KINDS_TYPES
function rvs_norm_0_${t1[0]}$${k1}$( ) result(res)
impure function rvs_norm_0_${t1[0]}$${k1}$ () result(res)
!
! Standard normal random vairate (0,1)
! Standard normal random variate (0,1)
!
${t1}$ :: res
${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$/r
Expand Down Expand Up @@ -157,8 +146,6 @@ contains

#:endfor



#:for k1, t1 in REAL_KINDS_TYPES
impure elemental &
function rvs_norm_${t1[0]}$${k1}$ (loc, scale) result(res)
Expand All @@ -168,19 +155,19 @@ contains
${t1}$, intent(in) :: loc, scale
${t1}$ :: res

if(scale == 0._${k1}$) call error_stop("Error(rvs_norm): Normal" &
//" distribution scale parameter must be non-zero")
if (scale <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = rvs_norm_0_${t1[0]}$${k1}$ ()
res = res*scale + loc
end if

end function rvs_norm_${t1[0]}$${k1}$

#:endfor



#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental &
function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res)
impure elemental function rvs_norm_${t1[0]}$${k1}$ (loc, scale) result(res)
!
! Normally distributed complex. The real part and imaginary part are &
! independent of each other.
Expand All @@ -192,24 +179,27 @@ contains
tr = rvs_norm_r${k1}$ (loc%re, scale%re)
ti = rvs_norm_r${k1}$ (loc%im, scale%im)
res = cmplx(tr, ti, kind=${k1}$)

end function rvs_norm_${t1[0]}$${k1}$

#:endfor



#:for k1, t1 in REAL_KINDS_TYPES
function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
impure function rvs_norm_array_${t1[0]}$${k1}$ (loc, scale, array_size) result(res)
${t1}$, intent(in) :: loc, scale
integer, intent(in) :: array_size
${t1}$ :: res(array_size)
${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$/r
${t1}$ :: x, y, re
integer :: hz, iz, i

if(scale == 0._${k1}$) call error_stop("Error(rvs_norm_array): Normal" &
//"distribution scale parameter must be non-zero")
if (.not. zig_norm_initialized) call zigset

if (scale <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
return
end if

do i = 1, array_size
iz = 0
hz = dist_rand(1_int32)
Expand Down Expand Up @@ -249,10 +239,8 @@ contains

#:endfor



#:for k1, t1 in CMPLX_KINDS_TYPES
function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
impure function rvs_norm_array_${t1[0]}$${k1}$ (loc, scale, array_size) result(res)
${t1}$, intent(in) :: loc, scale
integer, intent(in) :: array_size
integer :: i
Expand All @@ -264,33 +252,33 @@ contains
ti = rvs_norm_r${k1}$ (loc%im, scale%im)
res(i) = cmplx(tr, ti, kind=${k1}$)
end do

end function rvs_norm_array_${t1[0]}$${k1}$

#:endfor



#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
!
! Normal distribution probability density function
!
${t1}$, intent(in) :: x, loc, scale
${t1}$ :: res
${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$*acos(-1.0_${k1}$))

if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal" &
//"distribution scale parameter must be non-zero")
if (scale <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = exp(-0.5_${k1}$*((x - loc)/scale)*(x - loc)/scale)/ &
(sqrt_2_Pi*scale)
end if

end function pdf_norm_${t1[0]}$${k1}$

#:endfor



#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
elemental function pdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real(${k1}$) :: res

Expand All @@ -300,28 +288,27 @@ contains

#:endfor



#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
elemental function cdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
!
! Normal distribution cumulative distribution function
!
${t1}$, intent(in) :: x, loc, scale
${t1}$ :: res
${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$)

if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal" &
//"distribution scale parameter must be non-zero")
if (scale <= 0._${k1}$) then
res = ieee_value(1._${k1}$, ieee_quiet_nan)
else
res = erfc(-(x - loc)/(scale*sqrt_2))/2.0_${k1}$
end if

end function cdf_norm_${t1[0]}$${k1}$

#:endfor



#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res)
elemental function cdf_norm_${t1[0]}$${k1}$ (x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real(${k1}$) :: res

Expand Down