!   G13NBF Example Program Text
!   Mark 28.7 Release. NAG Copyright 2022.

    Module g13nbfe_mod

!     G13NBF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: costfn, get_data
    Contains
      Subroutine costfn(ts,nr,r,c,y,iuser,ruser,info)
!       Cost function, C. This cost function is based on the likelihood of
!       the gamma distribution

!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: info
        Integer, Intent (In)           :: nr, ts
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: c(nr)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*), y(0:*)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: r(nr)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: dn, shape, si
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: log, real
!       .. Executable Statements ..

!       RUSER(1) holds the shape parameter (a) for the gamma distribution
        shape = ruser(1)

!       Test which way around TS and R are (only needs to be done once)
        If (ts<r(1)) Then
          Do i = 1, nr
            si = y(r(i)) - y(ts)
            dn = real(r(i)-ts,kind=nag_wp)
            c(i) = 2.0_nag_wp*dn*shape*(log(si)-log(dn*shape))
          End Do

        Else
          Do i = 1, nr
            si = y(ts) - y(r(i))
            dn = real(ts-r(i),kind=nag_wp)
            c(i) = 2.0_nag_wp*dn*shape*(log(si)-log(dn*shape))
          End Do
        End If

!       Set info nonzero to terminate execution for any reason
        info = 0
      End Subroutine costfn

      Subroutine get_data(nin,n,k,y,iuser,ruser)
!       Read in data that is specific to the cost function

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: k
        Integer, Intent (In)           :: n, nin
!       .. Array Arguments ..
        Real (Kind=nag_wp), Allocatable, Intent (Out) :: ruser(:), y(:)
        Integer, Allocatable, Intent (Out) :: iuser(:)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: shape
        Integer                        :: i
!       .. Executable Statements ..

!       Read in the series of interest
!       NB: we are starting Y allocation at 0 as we manipulate
!       the data in Y in a moment
        Allocate (y(0:n))
        Read (nin,*) y(1:n)

!       Read in the shape parameter for the Gamma distribution
        Read (nin,*) shape

!       Store the shape parameter in RUSER. IUSER is not used
        Allocate (ruser(1),iuser(0))
        ruser(1) = shape

!       The cost function is a function of the sum of Y, so for
!       efficiency we will calculate the cumulative sum
!       It should be noted that this may introduce some rounding issues
!       with very extreme data
        y(0) = 0.0_nag_wp
        Do i = 1, n
          y(i) = y(i-1) + y(i)
        End Do

!       The value of K is defined by the cost function being used
!       in this example a value of 0.0 is the required value
        k = 0.0_nag_wp

        Return
      End Subroutine get_data
    End Module g13nbfe_mod

    Program g13nbfe

!     .. Use Statements ..
      Use g13nbfe_mod, Only: costfn, get_data
      Use nag_library, Only: g13nbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: beta, k
      Integer                          :: i, ifail, minss, n, ntau
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: ruser(:), y(:)
      Integer, Allocatable             :: iuser(:), tau(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: repeat
!     .. Executable Statements ..
      Write (nout,*) 'G13NBF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)

!     Read in the problem size, penalty and minimum segment size
      Read (nin,*) n, beta, minss

!     Read in the rest of the data, that (may be) dependent on the cost
!     function
      Call get_data(nin,n,k,y,iuser,ruser)

!     Allocate output arrays
      Allocate (tau(n))

!     Call routine to detect change points
      ifail = 0
      Call g13nbf(n,beta,minss,k,costfn,ntau,tau,y,iuser,ruser,ifail)

!     Display the results
      Write (nout,99999) ' -- Change Points --'
      Write (nout,99999) ' Number     Position'
      Write (nout,99999) repeat('=',21)
      Do i = 1, ntau
        Write (nout,99998) i, tau(i)
      End Do

99999 Format (1X,A)
99998 Format (1X,I4,7X,I6)
    End Program g13nbfe