! user_gamma 2 REAL REAL .TRUE.
! user_normal 2 REAL REAL .TRUE.
! user_defined_points 1 INTEGER .FALSE.
! user_defined_points2 2 INTEGER INTEGER .FALSE.
! user_bernoulli 1 REAL .FALSE.

MODULE user_dist
  IMPLICIT NONE

CONTAINS

  SUBROUTINE user_normal_gen(mu, s, buffer)
    ! A wrapper subroutine for the G05FDF NAG normal generation.
    IMPLICIT NONE
    INTEGER :: param
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), DIMENSION(:), INTENT(OUT) :: buffer
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(IN) :: mu, s
    EXTERNAL G05FDF

    CALL G05FDF (mu, s, SIZE(buffer), buffer)

  END SUBROUTINE user_normal_gen

  SUBROUTINE user_gamma_gen(a, b, buffer)
    IMPLICIT NONE
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(IN) :: a, b
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), DIMENSION(:), INTENT(OUT) :: buffer
    INTEGER :: ifail
    EXTERNAL G05FFF

    ifail = 0
    CALL g05FFF(a, b, SIZE(buffer), buffer, ifail)

    IF (ifail /= 0) THEN
       WRITE (*, *) 'Error when trying to generate proposals of the gamma distribution, parameters: a = ', a, ' b = ', b 
       STOP
    END IF

  END SUBROUTINE user_gamma_gen

  SUBROUTINE user_normal_freq(mu, sigma2, point, freq)
    ! calculates the logarithm of point's frequency in distribution Norm(mu, sigma) 
    IMPLICIT NONE
    INTRINSIC SQRT
    INTRINSIC LOG
    INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(8,0)
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: mu
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: sigma2
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: sigma
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN) :: point
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(OUT) :: freq

       sigma = SQRT(sigma2)
       freq = -LOG(sigma) + user_z_freq((point-mu)/sigma)       

  END SUBROUTINE user_normal_freq

  FUNCTION user_z_freq(x) RESULT (result)
    ! returns the logarithm of x's frequency in distribution Norm(0,1)
    IMPLICIT NONE
    INTRINSIC SQRT
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: x
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: result


    result =  -0.91893853 - 0.5*x*x

  END FUNCTION user_z_freq

  SUBROUTINE user_gamma_freq(alpha, lambda, point, freq)
    ! calculates the logarithm of point's frequency in distribution Gamma(alpha, lambda)
    IMPLICIT NONE
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: alpha, lambda, point
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(OUT):: freq

    freq =  alpha * LOG(lambda) + (alpha-1) * LOG(point) &
     - lambda * point - user_gammaln (alpha)

  END SUBROUTINE user_gamma_freq

  FUNCTION user_gammaln(xx) RESULT (result)
    ! returns the value of xx's gamma function
    IMPLICIT NONE
    INTRINSIC LOG
    REAL (KIND=SELECTED_REAL_KIND(8,0)), INTENT(IN):: xx
    REAL (KIND=SELECTED_REAL_KIND(8,0)) :: x, tmp, ser, result
    REAL (KIND=SELECTED_REAL_KIND(8,0)), DIMENSION(6) :: cof
    INTEGER :: j

    cof(1) =  76.18009173
    cof(2) = -86.50532033
    cof(3) =  24.01409822
    cof(4) = -1.231739516
    cof(5) =  0.120858003e-2
    cof(6) = -0.536382e-5

    x = xx - 1.0
    tmp = x + 5.5
    tmp = tmp - (x + 0.5) * LOG (tmp)
    ser = 1.0

    DO j = 1, 6
       x = x + 1.0;
       ser = ser + cof(j) / x;
    END DO

    result = -tmp + log (2.50662827465 * ser)

  END FUNCTION user_gammaln

  SUBROUTINE user_bernoulli_freq(p, point, prob)
    IMPLICIT NONE
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(IN) :: p
    INTEGER, INTENT(IN) :: point
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(OUT) :: prob


    IF (p < 0 .AND. point == 0) THEN
        prob = 1
        RETURN
    END IF

    IF (p < 0 .AND. point == 1) THEN
        prob = 0
        RETURN
    END IF

    IF (p > 1 .AND. point == 0) THEN
        prob = 0
        RETURN
    END IF

    IF (p > 1 .AND. point == 1) THEN
        prob = 1
        RETURN
    END IF

    IF (point == 1) THEN 
       prob = MIN(p, 1.0)
    ELSE
       IF (point == 0) THEN
          prob = 1-p
       ELSE
          prob = 0
       END IF
    END IF

  END SUBROUTINE user_bernoulli_freq

  SUBROUTINE user_defined_points_freq(param, point, prob)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: param
    INTEGER, INTENT(IN) :: point
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(OUT) :: prob

    ! P(obs = 0 | x = 0) = 0.8
    ! P(obs = 1 | x = 0) = 0.2
    ! P(obs = 0 | x = 1) = 0.4
    ! P(obs = 1 | x = 1) = 0.6

    IF (point == 0 .AND. param == 0) prob = 0.8
    IF (point == 1 .AND. param == 0) prob = 0.2
    IF (point == 0 .AND. param == 1) prob = 0.4
    IF (point == 1 .AND. param == 1) prob = 0.6

  END SUBROUTINE user_defined_points_freq

  SUBROUTINE user_defined_points2_freq(x, r, point, prob)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: x, r
    INTEGER, INTENT(IN) :: point
    REAL (KIND=SELECTED_REAL_KIND(8, 0)), INTENT(OUT) :: prob

    SELECT CASE (point)
    CASE (0)
       IF (x == 0) THEN
          ! P(obs=0 | x=0) = 0.8
          prob = 0.8
       ELSE
          SELECT CASE (r)
          CASE (1)
             ! P(obs=0 | x=1, r=1) = 0.3
             prob = 0.3
          CASE (2)
             ! P(obs=0 | x=1, r=2) = 0.2
             prob = 0.2
          CASE (3)
             ! P(obs=0 | x=1, r=3) = 0.1
             prob = 0.1
          CASE (4)
             ! P(obs=0 | x=1, r=4) = 0.05
             prob = 0.05
          CASE (5)
             ! P(obs=0 | x=1, r=5) = 0.05
             prob = 0.05
          END SELECT
       END IF
    CASE (1)
       IF (x == 0) THEN
          ! P(obs=1 | x=0) = 0.15
          prob = 0.15
       ELSE
          SELECT CASE (r)
          CASE (1)
             ! P(obs=1 | x=1, r=1) = 0.4
             prob = 0.4
          CASE (2)
             ! P(obs=1 | x=1, r=2) = 0.25
             prob = 0.25
          CASE (3)
             ! P(obs=1 | x=1, r=3) = 0.25
             prob = 0.25
          CASE (4)
             ! P(obs=1 | x=1, r=4) = 0.2
             prob = 0.2
          CASE (5)
             ! P(obs=1 | x=1, r=5) = 0.1
             prob = 0.1
          END SELECT

       END IF
    CASE (2)
       IF (x == 0) THEN
          ! P(obs=2 | x=0) = 0.05
          prob = 0.05
       ELSE
          SELECT CASE (r)
          CASE (1)
             ! P(obs=2 | x=1, r=1) = 0.2
             prob = 0.2
          CASE (2)
             ! P(obs=2 | x=1, r=2) = 0.3
             prob = 0.3
          CASE (3)
             ! P(obs=2 | x=1, r=3) = 0.3
             prob = 0.3
          CASE (4)
             ! P(obs=2 | x=1, r=4) = 0.25
             prob = 0.25
          CASE (5)
             ! P(obs=2 | x=1, r=5) = 0.25
             prob = 0.25
          END SELECT
       END IF

    CASE (3)
       IF (x == 0) THEN
          ! P(obs=3 | x=0) = 0.0
          prob = 0.0
       ELSE
          SELECT CASE (r)
          CASE (1)
             ! P(obs=3 | x=1, r=1) = 0.1
             prob = 0.1
          CASE (2)
             ! P(obs=3 | x=1, r=2) = 0.25
             prob = 0.25
          CASE (3)
             ! P(obs=3 | x=1, r=3) = 0.35
             prob = 0.35
          CASE (4)
             ! P(obs=3 | x=1, r=4) = 0.5
             prob = 0.5
          CASE (5)
             ! P(obs=3 | x=1, r=5) = 0.6
             prob = 0.6
          END SELECT ! select by r
       END IF ! choice by x
    END SELECT ! select by point

  END SUBROUTINE user_defined_points2_freq

END MODULE user_dist
