!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! functions.f90
!
! You are welcome to use this code for your own research but please 
! properly cite and attibute:
!
! Author: 
! Cameron Bracken, cameron.bracken@gmail.com
!
! Citation: 
! C. Bracken, K. Holman, B. Rajagopalan, H. Moradkhani, A Bayesian 
! hierarchical approach to multivariate nonstationary hydrologic frequency 
! analysis, Water Resources Research, 2017, Under Review. 
!
! Log probability density functions (and a few others) used in the 
! likelihood function. 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module functions
    implicit none

interface gev
  module procedure gev_s, gev_v
end interface

interface norm
  module procedure norm_s, norm_v, norm_m, norm_m3
end interface

interface uniform
  module procedure uniform_s, uniform_v, uniform_m
end interface

interface inv_logit
  module procedure inv_logit_s, inv_logit_v, inv_logit_m
end interface

interface std_norm
  module procedure std_norm_s, std_norm_v, std_norm_m, std_norm_m3
end interface

interface lognorm
  module procedure lognorm_s
end interface

    contains

pure function gev_s(y, loc, scale, shape)
    double precision:: gev_s
    double precision, intent(in)::y, loc, scale, shape
    double precision:: inv_shape, z

    inv_shape = 1.0/shape
    z = 1. + (y - loc) * shape / scale
    gev_s = -log(scale) - (1. + inv_shape)*log(z) - z**(-inv_shape)  
end function


pure function gev_v(y, loc, scale, shape)
    double precision:: gev_v
    double precision, intent(in):: y(:)
    double precision, intent(in):: loc, scale, shape
    integer:: n
    double precision, dimension(size(y)):: inv_shape, z

    n = size(y)

    inv_shape = 1.0/shape
    z = 1. + (y - loc) * shape / scale
    gev_v = -sum(log(scale) + (1. + inv_shape)*log(z) + z**(-inv_shape))
end function

pure function gev_cdf (q2, loc, scale, shape) 
    double precision, intent(in):: q2, loc, scale, shape 
    double precision:: q, p, r, gev_cdf

    q = (q2 - loc)/scale
    if (shape .eq. 0) then
        gev_cdf = exp(-exp(-q))
    else
        r = 1 + shape * q
        if(r < 0) r = 0.0
        gev_cdf = exp(-r**(-1/shape))
    end if

end function

function norm_s(x, mu, sigma)
    double precision, intent(in):: x, mu, sigma
    double precision:: norm_s
    !double precision, parameter:: pi=3.141592653589793238462643d0


    !norm_s = -log(sigma) - 0.5*log(2*pi) - 0.5*(x-mu)*(x-mu) / (sigma*sigma)
    ! dropping constant term
    norm_s = -log(sigma) - 0.5*(x-mu)*(x-mu) / (sigma*sigma)
end function

function norm_v(x, mu, sigma)
    double precision, intent(in):: x(:)
    double precision, intent(in):: mu, sigma
    double precision:: norm_v
    !double precision, parameter:: pi=3.141592653589793238462643d0
    integer:: n

    n = size(x)

    !norm_v = -n *( log(sigma) + 0.5*log(2.0*pi) ) - &
    !        0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
    ! dropping constant term
    norm_v = -n * log(sigma) - 0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
end function

function norm_m(x, mu, sigma)
    double precision, intent(in):: x(:,:)
    double precision, intent(in):: mu, sigma
    double precision:: norm_m
    !double precision, parameter:: pi=3.141592653589793238462643d0
    integer:: n1, n2

    n1 = size(x,1)
    n2 = size(x,1)

    !norm_v = -n *( log(sigma) + 0.5*log(2.0*pi) ) - &
    !        0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
    ! dropping constant term
    norm_m = -n1 * n2 * log(sigma) - 0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
end function

function norm_m3(x, mu, sigma)
    double precision, intent(in):: x(:,:,:)
    double precision, intent(in):: mu, sigma
    double precision:: norm_m3
    !double precision, parameter:: pi=3.141592653589793238462643d0
    integer:: n1, n2

    n1 = size(x,1)
    n2 = size(x,1)

    !norm_v = -n *( log(sigma) + 0.5*log(2.0*pi) ) - &
    !        0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
    ! dropping constant term
    norm_m3 = -n1 * n2 * log(sigma) - 0.5 * sum((x-mu)*(x-mu)) / (sigma*sigma)
end function


function lognorm_s(x, mu, sigma)
    double precision, intent(in):: x, mu, sigma
    double precision:: lognorm_s, lx
    !double precision, parameter:: pi=3.141592653589793238462643d0


    lx = log(x)
    ! dropping constant term
    lognorm_s = - lx - log(sigma) - 0.5*(lx-mu)*(lx-mu) / (sigma*sigma)
end function

function uniform_s(x, a, b)
    double precision, intent(in):: x, a, b
    double precision:: uniform_s
    double precision, parameter:: inf = 1.7976931348623157d308

    uniform_s = -inf
    if(x < a .or. x > b) return
    uniform_s = 1./ (b-a)
end function

function uniform_v(x, a, b)
    double precision, intent(in):: x(:), a, b
    double precision:: uniform_v
    double precision, parameter:: inf = 1.7976931348623157d308
    integer:: i

    uniform_v = -inf
    if(any(x < a) .or. any(x > b)) return
    uniform_v = size(x) / (b-a)
end function

function uniform_m(x, a, b)
    double precision, intent(in):: x(:,:), a, b
    double precision:: uniform_m
    double precision, parameter:: inf = 1.7976931348623157d308
    integer:: i,j

    uniform_m = -inf
    if(any(x < a) .or. any(x > b)) return
    uniform_m = size(x,1) * size(x,2) / (b-a)
end function

!function norm_s(x, mu, sigma)
!    integer:: n, nmu, nsigma
!    double precision:: like, norm_s
!    double precision, intent(in):: x, mu, sigma
!
!    nsigma = 1
!    nmu = 1
!    n = 1
!    call normal(x, mu, sigma, n, nmu, nsigma, like)
!
!    norm_s = like
!    write(*,*)like
!end function
!
!function norm_v(x, mu, sigma)
!    integer:: nmu, nsigma
!    double precision:: like, norm_v
!    double precision, intent(in):: x(:)
!    double precision, intent(in):: mu, sigma
!    integer:: n
!
!    n = size(x)
!
!    nsigma = 1
!    nmu = 1
!    call normal(x, mu, sigma, n, nmu, nsigma, like)
!
!    norm_v = like
!    write(*,*)like
!end function

function std_norm_s(x)
    integer:: n
    double precision:: std_norm_s
    double precision, intent(in):: x
    double precision:: mu, sigma

    mu = 0.
    sigma = 1.
    std_norm_s = norm(x, mu, sigma)

    !std_norm_s = like
    !write(*,*)like
end function

function std_norm_v(x)
    double precision:: like, std_norm_v
    double precision:: x(:)
    double precision:: mu, sigma


    mu = 0.
    sigma = 1.
    std_norm_v = norm(x, mu, sigma)

end function

function std_norm_m(x)
    double precision:: like, std_norm_m
    double precision:: x(:,:)
    double precision:: mu, sigma


    mu = 0.
    sigma = 1.
    std_norm_m = norm(x, mu, sigma)

end function

function std_norm_m3(x)
    double precision:: like, std_norm_m3
    double precision:: x(:,:,:)
    double precision:: mu, sigma


    mu = 0.
    sigma = 1.
    std_norm_m3 = norm(x, mu, sigma)

end function

function inv_logit_s(x)
    double precision, intent(in):: x
    double precision:: inv_logit_s
    inv_logit_s = 1./(1.+exp(-x))
end function

function inv_logit_v(x)
    double precision, dimension(:), intent(in):: x
    double precision, dimension(size(x)):: inv_logit_v
    inv_logit_v = 1./(1.+exp(-x))
end function

function inv_logit_m(x)
    double precision, intent(in):: x(:,:)
    double precision:: inv_logit_m(size(x,1),size(x,2))
    inv_logit_m = 1./(1.+exp(-x))
end function

function isinf(x)
    logical:: isinf
    double precision, intent(in)::x
    double precision, parameter:: infinity = 1.7976931348623157d308
    if(abs(x) > infinity) then
        isinf = .true.
    else
        isinf = .false.
    end if
end function

function exp_cov_matrix (x, sill, range, nugget)

    double precision, intent(in):: x(:,:), sill, range, nugget
    double precision:: exp_cov_matrix(size(x,1), size(x,1))
    double precision:: neg_inv_range
    integer:: i,j, S

    S = size(x,1)
    neg_inv_range = -1.d0/range

    !write(*,*) sill, range, neg_inv_range

    ! off-diagonal elements
    do i = 1,(S-1)
      do j = (i+1),S
        !write(*,*)sill * exp(neg_inv_range * x(i,j))
        exp_cov_matrix(i,j) = sill * exp(neg_inv_range * x(i,j))
        exp_cov_matrix(j,i) = exp_cov_matrix(i,j)
      end do
    end do

    ! diagonal elements
    do i = 1,S
      exp_cov_matrix(i,i) = sill + nugget
    end do

end function



function gauss_eliptical_copula(logf, F, D)

    double precision:: logf(:), F(:), D(:,:)
    double precision, dimension(size(F)):: u, zeros
    integer i, n, info
    double precision:: like, gauss_eliptical_copula, sum_log_phi_inv_u

    n = size(F)
    zeros = 0.

    !F[F==1] = 0.999999
    do i = 1,n
        if(F(i) == 1)F(i) = 0.99999
        if(F(i) == 0)F(i) = 0.00001
        u(i) = inv_norm_cdf(F(i))
        !u = r8_normal_01_cdf_inverse(F(i))
    end do
    !write(*,*) F(1),u(1)

    sum_log_phi_inv_u = std_norm(u)

    call chol_mvnorm(u, zeros, D, n, like, info)

    gauss_eliptical_copula = sum(logf) - sum_log_phi_inv_u + like

    !write(*,*)sum(logf),sum_log_phi_inv_u
    !write(*,*)F
    !write(*,*)logf
    !write(*,*)sum(logf), sum(log_phi_u), like
    
end function

function inv_norm_cdf(p)
    double precision, intent(in):: p
    double precision:: x, q, r, p_low, p_high, e, u, inv_norm_cdf
    double precision:: a(6), b(5), c(6), d(4)
    double precision, parameter :: pi = 4 * atan (1d0)

    a(1) = -3.969683028665376e+01
    a(2) =  2.209460984245205e+02
    a(3) = -2.759285104469687e+02
    a(4) =  1.383577518672690e+02
    a(5) = -3.066479806614716e+01
    a(6) =  2.506628277459239e+00

    b(1) = -5.447609879822406e+01
    b(2) =  1.615858368580409e+02
    b(3) = -1.556989798598866e+02
    b(4) =  6.680131188771972e+01
    b(5) = -1.328068155288572e+01

    c(1) = -7.784894002430293e-03
    c(2) = -3.223964580411365e-01
    c(3) = -2.400758277161838e+00
    c(4) = -2.549732539343734e+00
    c(5) =  4.374664141464968e+00
    c(6) =  2.938163982698783e+00

    d(1) =  7.784695709041462e-03
    d(2) =  3.224671290700398e-01
    d(3) =  2.445134137142996e+00
    d(4) =  3.754408661907416e+00

    if ( p <= 0. ) then
        inv_norm_cdf = - huge ( p )
        return
    end if

    if ( p >= 1. ) then
        inv_norm_cdf = huge ( p )
        return
    end if

    p_low  = 0.02425
    p_high = 0.97575 ! 1 - p_low
    if ( (p_low <= p) .and. (p <= p_high) ) then ! central region
      q = p - 0.5
      r = q * q
      x = (((((a(1)*r+a(2))*r+a(3))*r+a(4))*r+a(5))*r+a(6))*q / &
           (((((b(1)*r+b(2))*r+b(3))*r+b(4))*r+b(5))*r+1.0)      
    else if (p < p_low) then ! lower region
      q = sqrt(-2.0*log(p))
      x = (((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /  &
            ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1.0)            
    else  ! upper region
      q = sqrt(-2.0*log(1-p))
      x = -(((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /  &
             ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1.0)
    end if

    ! refinement
    e = norm_cdf(x) - p
    u = e * sqrt(2.0 * pi) * exp(0.5 * x * x)
    x = x - u / (1.0 + 0.5 * x * u)
    inv_norm_cdf = x
end function


function norm_cdf(x)
    double precision:: norm_cdf, x, cdf, ccdf

    call norm_cdf_sub(x, cdf, ccdf)
    norm_cdf = cdf

end function

function lognorm_cdf(x, mu, sigma)
    double precision:: lognorm_cdf, x, mu, sigma

    lognorm_cdf = norm_cdf((log(x) - mu)/sigma)

end function


function mean(x)
  double precision, intent(in):: x(:)
  double precision:: mean
  integer:: n
  
  n = size(x)

  mean = sum(x)/dble(n)
end function

function sd(x)

  double precision, intent(in):: x(:)
  double precision:: sd, xbar
  integer:: n
  
  n = size(x)

  xbar = mean(x)
  sd = sqrt(sum((x-xbar)**2d0)/dble(n-1))

end function 

! function rhat(ch)
!   double precision:: rhat, ch(:,:), W, B, nd
!   double precision:: sds(size(ch,2)), sds(size(ch,2)), chain(size(ch,1)), 
!   integer:: m, n, i

!   n = size(ch,1)
!   nd = dble(n)
!   m = size(ch,2)

!   do i = 1,m
!     do j = 1,n
!       chain(j) = ch(j,i)
!     end do 
!     sds(i) = sd(chain)
!     means(i) = mean(chain)
!   end do 

!   W = mean(sds)
!   B = sd(means)**2d0

!   rhat = sqrt(((1-1/nd) * W + 1/nd * B)/W)

! end function 

subroutine norm_cdf_sub ( arg, cum, ccum )

!*****************************************************************************80
!
!! CUMNOR computes the cumulative normal distribution.
!
!  Discussion:
!
!    This function evaluates the normal distribution function:
!
!                              / x
!                     1       |       -t*t/2
!          P(x) = ----------- |      e       dt
!                 sqrt(2 pi)  |
!                             /-oo
!
!    This transportable program uses rational functions that
!    theoretically approximate the normal distribution function to
!    at least 18 significant decimal digits.  The accuracy achieved
!    depends on the arithmetic system, the compiler, the intrinsic
!    functions, and proper selection of the machine dependent
!    constants.
!
!  Author:
!
!    William Cody
!    Mathematics and Computer Science Division
!    Argonne National Laboratory
!    Argonne, IL 60439
!
!  Reference:
!
!    William Cody,
!    Rational Chebyshev approximations for the error function,
!    Mathematics of Computation,
!    1969, pages 631-637.
!
!    William Cody,
!    Algorithm 715:
!    SPECFUN - A Portable FORTRAN Package of Special Function Routines
!    and Test Drivers,
!    ACM Transactions on Mathematical Software,
!    Volume 19, Number 1, 1993, pages 22-32.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ARG, the upper limit of integration.
!
!    Output, real ( kind = 8 ) CUM, CCUM, the Normal density CDF and
!    complementary CDF.
!
!  Local Parameters:
!
!    Local, real ( kind = 8 ) EPS, the argument below which anorm(x)
!    may be represented by 0.5 and above which  x*x  will not underflow.
!    A conservative value is the largest machine number X
!    such that   1.0D+00 + X = 1.0D+00   to machine precision.
!
  implicit none

  real ( kind = 8 ), parameter, dimension ( 5 ) :: a = (/ &
    2.2352520354606839287D+00, &
    1.6102823106855587881D+02, &
    1.0676894854603709582D+03, &
    1.8154981253343561249D+04, &
    6.5682337918207449113D-02 /)
  real ( kind = 8 ) arg
  real ( kind = 8 ), parameter, dimension ( 4 ) :: b = (/ &
    4.7202581904688241870D+01, &
    9.7609855173777669322D+02, &
    1.0260932208618978205D+04, &
    4.5507789335026729956D+04 /)
  real ( kind = 8 ), parameter, dimension ( 9 ) :: c = (/ &
    3.9894151208813466764D-01, &
    8.8831497943883759412D+00, &
    9.3506656132177855979D+01, &
    5.9727027639480026226D+02, &
    2.4945375852903726711D+03, &
    6.8481904505362823326D+03, &
    1.1602651437647350124D+04, &
    9.8427148383839780218D+03, &
    1.0765576773720192317D-08 /)
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ), parameter, dimension ( 8 ) :: d = (/ &
    2.2266688044328115691D+01, &
    2.3538790178262499861D+02, &
    1.5193775994075548050D+03, &
    6.4855582982667607550D+03, &
    1.8615571640885098091D+04, &
    3.4900952721145977266D+04, &
    3.8912003286093271411D+04, &
    1.9685429676859990727D+04 /)
  real ( kind = 8 ) del
  real ( kind = 8 ) eps
  integer ( kind = 4 ) i
  real ( kind = 8 ), parameter, dimension ( 6 ) :: p = (/ &
    2.1589853405795699D-01, &
    1.274011611602473639D-01, &
    2.2235277870649807D-02, &
    1.421619193227893466D-03, &
    2.9112874951168792D-05, &
    2.307344176494017303D-02 /)
  real ( kind = 8 ), parameter, dimension ( 5 ) :: q = (/ &
    1.28426009614491121D+00, &
    4.68238212480865118D-01, &
    6.59881378689285515D-02, &
    3.78239633202758244D-03, &
    7.29751555083966205D-05 /)
  real ( kind = 8 ), parameter :: root32 = 5.656854248D+00
  real ( kind = 8 ), parameter :: sixten = 16.0D+00
  real ( kind = 8 ) temp
  real ( kind = 8 ), parameter :: sqrpi = 3.9894228040143267794D-01
  real ( kind = 8 ), parameter :: thrsh = 0.66291D+00
  real ( kind = 8 ) x
  real ( kind = 8 ) xden
  real ( kind = 8 ) xnum
  real ( kind = 8 ) y
  real ( kind = 8 ) xsq
!
!  Machine dependent constants
!
  eps = epsilon ( 1.0D+00 ) * 0.5D+00

  x = arg
  y = abs ( x )

  if ( y <= thrsh ) then
!
!  Evaluate  anorm  for  |X| <= 0.66291
!
    if ( eps < y ) then
      xsq = x * x
    else
      xsq = 0.0D+00
    end if

    xnum = a(5) * xsq
    xden = xsq
    do i = 1, 3
      xnum = ( xnum + a(i) ) * xsq
      xden = ( xden + b(i) ) * xsq
    end do
    cum = x * ( xnum + a(4) ) / ( xden + b(4) )
    temp = cum
    cum = 0.5D+00 + temp
    ccum = 0.5D+00 - temp
!
!  Evaluate ANORM for 0.66291 <= |X| <= sqrt(32)
!
  else if ( y <= root32 ) then

    xnum = c(9) * y
    xden = y
    do i = 1, 7
      xnum = ( xnum + c(i) ) * y
      xden = ( xden + d(i) ) * y
    end do
    cum = ( xnum + c(8) ) / ( xden + d(8) )
    xsq = aint ( y * sixten ) / sixten
    del = ( y - xsq ) * ( y + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) * exp ( -del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if
!
!  Evaluate ANORM for sqrt(32) < |X|.
!
  else

    cum = 0.0D+00
    xsq = 1.0D+00 / ( x * x )
    xnum = p(6) * xsq
    xden = xsq
    do i = 1, 4
      xnum = ( xnum + p(i) ) * xsq
      xden = ( xden + q(i) ) * xsq
    end do

    cum = xsq * ( xnum + p(5) ) / ( xden + q(5) )
    cum = ( sqrpi - cum ) / y
    xsq = aint ( x * sixten ) / sixten
    del = ( x - xsq ) * ( x + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) &
      * exp ( - del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if

  end if

  if ( cum < tiny ( cum ) ) then
    cum = 0.0D+00
  end if

  if ( ccum < tiny ( ccum ) ) then
    ccum = 0.0D+00
  end if

  return
end subroutine

subroutine r8_swap ( x, y )
  real ( kind = 8 ) x, y, z

  z = x
  x = y
  y = z

  return
end subroutine 

subroutine init_random_seed()
  use iso_fortran_env, only: int64
  implicit none
  integer, allocatable :: seed(:)
  integer :: i, n, un, istat, dt(8), pid, getpid
  integer(int64) :: t

  call random_seed(size = n)
  allocate(seed(n))
  ! First try if the OS provides a random number generator
  open(newunit=un, file="/dev/urandom", access="stream", &
       form="unformatted", action="read", status="old", iostat=istat)
  if (istat == 0) then
     read(un) seed
     close(un)
  else
     ! Fallback to XOR:ing the current time and pid. The PID is
     ! useful in case one launches multiple instances of the same
     ! program in parallel.
     call system_clock(t)
     if (t == 0) then
        call date_and_time(values=dt)
        t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
             + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
             + dt(3) * 24_int64 * 60 * 60 * 1000 &
             + dt(5) * 60 * 60 * 1000 &
             + dt(6) * 60 * 1000 + dt(7) * 1000 &
             + dt(8)
     end if
     pid = getpid()
     t = ieor(t, int(pid, kind(t)))
     do i = 1, n
        seed(i) = lcg(t)
     end do
  end if
  call random_seed(put=seed)
contains
  ! This simple PRNG might not be good enough for real work, but is
  ! sufficient for seeding a better PRNG.
  function lcg(s)
    integer :: lcg
    integer(int64) :: s
    if (s == 0) then
       s = 104729
    else
       s = mod(s, 4294967296_int64)
    end if
    s = mod(s * 279470273_int64, 4294967291_int64)
    lcg = int(mod(s, int(huge(0), int64)), kind(0))
  end function lcg
end subroutine init_random_seed


end module