!   H05ABF Example Program Text

!   Mark 24 Release. NAG Copyright 2012.
    Module h05abfe_mod

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: nin = 5, nout = 6
    Contains
      Subroutine f(m,drop,lz,z,la,a,score,iuser,ruser,info)
!       Score calculating function required by H05ABF

!       M,DROP,LZ,Z,LA,A,SCORE IUSER and RUSER are all as described in the 
!       documentation of H05ABF.

!       This particular example finds the set, of a given size, of explanatory 
!       variables that best fit a response variable when a linear regression 
!       model is used. Therefore the feature set is the set of all the 
!       explanatory variables and the best set of features is defined as set of
!       explanatory variables that gives the smallest residual sums of squares.
!       See the documentation for G02DAF for details on linear regression
!       models.

!       .. Use Statements ..
        Use nag_library, Only: g02daf
!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: drop, la, lz, m
        Integer, Intent (Inout)              :: info
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: ruser(*)
        Real (Kind=nag_wp), Intent (Out)     :: score(max(la,1))
        Integer, Intent (In)                 :: a(la), z(lz)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: rss, tol
        Integer                              :: ex, ey, i, idf, ifail,         &
                                                inv_drop, ip, irank, ldq, ldx, &
                                                n, sx, sy
        Logical                              :: svd
        Character (1)                        :: mean, weight
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable      :: b(:), cov(:), h(:), p(:),      &
                                                q(:,:), res(:), se(:), wk(:),  &
                                                wt(:)
        Integer, Allocatable                 :: isx(:)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs, count, max
!       .. Executable Statements ..
        Continue
        info = 0

        n = iuser(1)
        ldq = n
        ldx = n

!       No intercept term and no weights
        mean = 'Z'
        weight = 'U'

!       Allocate various arrays required by G02DAF
        Allocate (b(m),cov((m*m+m)/2),h(n),p(m*(m+2)),q(ldq,m+1),res(n),se(m), &
          wk(m*m+5*(m-1)),wt(1),isx(m))

!       Set up the initial feature set. If DROP = 0, this is the Null set (i.e.
!       no features). If DROP = 1 then this is the full set (i.e. all features)
        isx(1:m) = drop

!       Add (if DROP = 0) or remove (if DROP = 1) the all the features specified
!       in Z
        inv_drop = abs(drop-1)
        Do i = 1, lz
          isx(z(i)) = inv_drop
        End Do

!       Get the start and end of X and Y in RUSER
        sx = iuser(4)
        ex = sx + m*n - 1
        sy = iuser(3)
        ey = sy + n - 1

!       Extract some parameters from RUSER
        tol = ruser(iuser(2))

        Do i = 1, max(la,1)
          If (la>0) Then
            If (i>1) Then
!             Reset the feature altered at the last iteration
              isx(a(i-1)) = drop
            End If

!           Add or drop the I'th feature in A 
            isx(a(i)) = inv_drop
          End If

          ip = count(isx(1:m)==1)

!         Fit the regression model
          ifail = 0
          Call g02daf(mean,weight,n,ruser(sx:ex),ldx,m,isx,ip,ruser(sy:ey),wt, &
            rss,idf,b,se,cov,res,h,q,ldq,svd,irank,p,tol,wk,ifail)

!         Return the score (the residual sums of squares)
          score(i) = rss
        End Do

!       Keep track of the number of subsets evaluated
        iuser(5) = iuser(5) + max(1,la)
      End Subroutine f

      Subroutine prepare_user_arrays(m,iuser,ruser)
!       Populate the user arrays

!       M is the maximum number of features (as per H05ABF)

!       In this example RUSER holds the data required for the linear regression
!       (the matrix of explanatory variables, X, the vector of response values,
!       Y, and the tolerance, TOL) and IUSER holds number of observations and 
!       the location of the various elements in RUSER

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: m
!       .. Array Arguments ..
        Real (Kind=nag_wp), Allocatable, Intent (Out) :: ruser(:)
        Integer, Allocatable, Intent (Out)   :: iuser(:)
!       .. Local Scalars ..
        Integer                              :: i, ierr, j, liuser, lruser, n, &
                                                p1, p2
        Character (200)                      :: line
!       .. Executable Statements ..
        Continue

!       Read in the number of observations for the data used in the linear 
!       regression skipping any headings or blank lines
        Do
          Read (nin,*,Iostat=ierr) line
          If (ierr/=0) Exit
          Read (line,*,Iostat=ierr) n
          If (ierr==0) Exit
        End Do

!       Allocate space for the user arrays
        liuser = 5
        lruser = n + n*m + 1
        Allocate (ruser(lruser),iuser(liuser))

!       Number of observations
        iuser(1) = n
!       Location of TOL in RUSER
        iuser(2) = 1
!       Start of Y in RUSER
        iuser(3) = 2
!       Start of X in RUSER
        iuser(4) = iuser(3) + n
!       Keep track of the number of subsets evaluated
        iuser(5) = 0

!       Read in the tolerance for the regression
        Read (nin,*) ruser(iuser(2))

!       Read in the data
        p1 = iuser(3)
        p2 = iuser(4)
        Do i = 0, n - 1
          Read (nin,*)(ruser(p2+i+j*n),j=0,m-1), ruser(p1+i)
        End Do
      End Subroutine prepare_user_arrays
    End Module h05abfe_mod

    Program h05abfe

!     .. Use Statements ..
      Use nag_library, Only: h05abf, nag_wp
      Use h05abfe_mod, Only: f, nin, nout, prepare_user_arrays
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: gamma
      Integer                              :: i, ifail, ip, j, la, m, mincnt,  &
                                              mincr, mip, nbest
!     .. Local Arrays ..
      Real (Kind=nag_wp)                   :: acc(2)
      Real (Kind=nag_wp), Allocatable      :: bscore(:), ruser(:)
      Integer, Allocatable                 :: a(:), bz(:,:), ibz(:), iuser(:), &
                                              z(:)
      Logical, Allocatable                 :: mask(:)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max, pack
!     .. Executable Statements ..
      Write (nout,*) 'H05ABF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) m, ip, nbest

!     Read in the control parameters for the subset selection
      Read (nin,*) mincr, mincnt, gamma, acc(1:2)

!     Allocate memory required by the subset selection routine
      mip = m - ip
      Allocate (z(mip),a(max(nbest,m)),bz(mip,nbest),bscore(max(nbest,m)))

!     Prepare the user workspace arrays IUSER and RUSER
      Call prepare_user_arrays(m,iuser,ruser)

!     Call the forward communication best subset routine
      ifail = -1
      Call h05abf(mincr,m,ip,nbest,la,bscore,bz,f,mincnt,gamma,acc,iuser, &
        ruser,ifail)
      If (ifail/=0 .And. ifail/=42) Then
!       An error occurred
        Go To 100
      End If

!     Titles
      Write (nout,99999) '    Score        Feature Subset'
      Write (nout,99999) '    -----        --------------'

!     Display the best subsets and corresponding scores. H05AAF returns a list 
!     of features excluded from the best subsets, so this is inverted to give
!     the set of features included in each subset
      Allocate (ibz(m),mask(m))
      ibz(1:m) = (/(i,i=1,m)/)
      Do i = 1, la
        mask(1:m) = .True.
        Do j = 1, mip
          mask(bz(j,i)) = .False.
        End Do
        Write (nout,99998) bscore(i), pack(ibz,mask)
      End Do

      Write (nout,*)
      If (ifail==42) Then
        Write (nout,99997) nbest, &
          ' subsets of the requested size do not exist, only ', la, &
          ' are displayed.'
      End If
      Write (nout,99996) iuser(5), ' subsets evaluated in total'

100   Continue

99999 Format (1X,A)
99998 Format (1X,E12.5,100(1X,I5))
99997 Format (1X,I0,A,I0,A)
99996 Format (1X,I0,A)

    End Program h05abfe
