NAG Library Manual, Mark 28.5
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   G13NEF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.
    Module g13nefe_mod

!     G13NEF 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                           :: chgpfn, get_data
    Contains
      Subroutine chgpfn(side,u,w,minss,v,cost,y,iuser,ruser,info)
!       Routine to calculate a proposed change point and associated cost
!       The cost is based on the likelihood of the gamma distribution

!       .. Use Statements ..
        Use nag_library, Only: x07caf, x07cbf
!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: info
        Integer, Intent (In)           :: minss, side, u, w
        Integer, Intent (Out)          :: v
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: cost(3)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(0:*), y(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: dn, shape, this_cost, tmp, ys
        Integer                        :: floc, i, li, lloc
!       .. Local Arrays ..
        Integer                        :: cexmode(3), texmode(3)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: log
!       .. Executable Statements ..

!       The gamma cost function used below can result in log(0) being taken
!       (if there is a segment of zeros in Y), this leads to a cost of -Inf
!       (which is correct), but we need to make sure that the compiler
!       doesn't stop at the creation of the -Inf

!       Save the current IEEE exception mode
        Call x07caf(cexmode)

!       Set the IEEE exception mode to not trap division by zero
        texmode(:) = cexmode(:)
        texmode(2) = 0
        Call x07cbf(texmode)

!       Extract shape from RUSER
        shape = ruser(0)

!       Calculate the first and last positions for potential change
!       points, conditional on the fact that each sub-segment must be
!       at least MINSS wide
        floc = u + minss - 1
        lloc = w - minss

!       In order to calculate the cost of having a change point at I, we
!       need to calculate C(Y(FLOC:I)) and C(Y(I+1:LLOC)), where C(.) is
!       the cost function (based on the gamma distribution in this example).
!       Rather than calculate these values at each call to CHGPFN we store
!       the values for later use

!       If SIDE = 1 (i.e. we are working with a left hand sub-segment),
!       we already have C(Y(FLOC:I)) for this value of FLOC, so only need
!       to calculate C(Y(I+1:LLOC)), similarly when SIDE = 2 we only need
!       to calculate C(Y(FLOC:I))
!       When SIDE = -1, we need the cost of the full segment, which we do
!       in a forwards manner (calculating C(Y(FLOC:I)) in the process), so
!       when SIDE = 0 we only need to calculate C(Y(I:1:LLOC))

!       Get the intermediate costs
        ys = 0.0_nag_wp
        dn = 0.0_nag_wp
        If (side==0 .Or. side==1) Then
!         RUSER(2*I) = C(Y(I+1:W))
          Do i = w, floc + 1, -1
            dn = dn + 1.0_nag_wp
            tmp = dn*shape
            ys = ys + y(i)
            ruser(2*i-2) = 2.0_nag_wp*tmp*(log(ys)-log(tmp))
          End Do

        Else
!         RUSER(2*I-1) = C(Y(U:I))
          If (side==-1) Then
            li = w
          Else
            li = lloc
          End If
          Do i = u, li
            dn = dn + 1.0_nag_wp
            tmp = dn*shape
            ys = ys + y(i)
            ruser(2*i-1) = 2.0_nag_wp*tmp*(log(ys)-log(tmp))
          End Do
        End If

        If (side>=0) Then
!         Need to find a potential change point
          v = 0
          cost(1) = 0.0_nag_wp

!         Loop over all possible change point locations
          Do i = floc, lloc
            this_cost = ruser(2*i-1) + ruser(2*i)

            If (this_cost<cost(1) .Or. v==0) Then
!             Update the proposed change point location
              v = i
              cost(1) = this_cost
              cost(2) = ruser(2*i-1)
              cost(3) = ruser(2*i)
            End If
          End Do
        Else
!         Need to calculate the cost for the full segment
          cost(1) = ruser(2*w-1)
!         No need to populate the rest of COST or V
        End If

!       Reset the IEEE exception mode back to what it was
        Call x07cbf(cexmode)

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

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

!       .. Scalar Arguments ..
        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
!       .. Executable Statements ..

!       Read in the series of interest
        Allocate (y(1:n))
        Read (nin,*) y(1:n)

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

!       We are going to use RUSER for two purposes, firstly to store the shape
!       parameter, and we also need an additional 2*N elements of workspace
!       we reference from 0 to make the coding easier later

!       IUSER is not going to be used
        Allocate (iuser(0),ruser(0:2*n))

!       Store the shape parameter in the 0th element of RUSER
        ruser(0) = shape

!       We will be populating the other elements of RUSER in the first
!       call to CHGPFN

        Return
      End Subroutine get_data
    End Module g13nefe_mod

    Program g13nefe

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

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

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

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

!     Allocate output arrays
      Allocate (tau(n))

!     Call routine to detect change points
      ifail = 0
      Call g13nef(n,beta,minss,mdepth,chgpfn,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 g13nefe