! G13NEF Example Program Text
! Mark 27.1 Release. NAG Copyright 2020.
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 ..
Continue
! 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 ..
Continue
! 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 ..
Continue
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