!   G22YAF Example Program Text
!   Mark 29.1 Release. NAG Copyright 2023.

    Module g22yafe_mod
!     G22YAF Example Program Module:
!            Parameters and User-defined Routines

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: read_line
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6

    Contains
      Subroutine read_line(v1)
!       Read in a line from NIN and remove any comments

!       .. Scalar Arguments ..
        Character (*), Intent (Out)    :: v1
!       .. Local Scalars ..
        Integer                        :: pend
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, index
!       .. Executable Statements ..

        Read (nin,'(A200)') v1
        pend = index(v1,'::')
        If (pend/=0) Then
          v1 = v1(1:pend-1)
        End If
        v1 = adjustl(v1)

        Return
      End Subroutine read_line
    End Module g22yafe_mod

    Program g22yafe

!     .. Use Statements ..
      Use g22yafe_mod, Only: nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g22yaf, g22ybf, g22ycf, g22zaf, g22zmf, g22znf,   &
                             nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hform, hxdesc
      Real (Kind=nag_wp)               :: rvalue
      Integer                          :: i, ifail, ivalue, lddat, ldx,        &
                                          lvnames, mx, nobs, nvar, optype,     &
                                          sddat, sdx
      Character (200)                  :: cvalue, formula
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dat(:,:), x(:,:), y(:)
      Real (Kind=nag_wp)               :: tx(0,0)
      Integer, Allocatable             :: levels(:)
      Character (50), Allocatable      :: vnames(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: trim
!     .. Executable Statements ..
      Write (nout,*) 'G22YAF Example Program Results'
      Write (nout,*)

      hform = c_null_ptr
      hddesc = c_null_ptr
      hxdesc = c_null_ptr

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

!     Read in the model formula, remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hform,formula,ifail)

!     Extract and display the parsed formula
      ifail = 0
      Call g22znf(hform,'Formula',ivalue,rvalue,cvalue,optype,ifail)
      Write (nout,*) 'Formula: ', trim(cvalue)
      Write (nout,*)

!     Read in size of the data matrix and number of variable labels supplied
      Read (nin,*) nobs, nvar, lvnames

!     Read in number of levels and names for the variables
      Allocate (levels(nvar),vnames(lvnames))
      Read (nin,*) levels(1:nvar)
      If (lvnames>0) Then
        Read (nin,*) vnames(1:lvnames)
      End If

!     Create a description of the data matrix
      ifail = 0
      Call g22ybf(hddesc,nobs,nvar,levels,lvnames,vnames,ifail)

!     Read in the data matrix and response variable
      lddat = nobs
      sddat = nvar
      Allocate (dat(lddat,sddat),y(nobs))
      Read (nin,*)(dat(i,1:nvar),y(i),i=1,nobs)

!     We want the design matrix to include an explicit term for the mean effect
      ifail = 0
      Call g22zmf(hform,'Explicit Mean = Yes',ifail)

!     Calculate the size of the design matrix
      ldx = 0
      sdx = 0
      ifail = 1
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      If (ifail/=81 .And. ifail/=82 .And. ifail/=91 .And. ifail/=92) Then
!       Redisplay any error messages not related to X being too small
        ifail = 0
        Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      End If

!     Generate the design matrix
      ldx = nobs
      sdx = mx
      Allocate (x(ldx,sdx))
      ifail = 0
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,x,ldx,sdx,mx,ifail)

!     Display the design matrix
      Write (nout,*) 'Design Matrix (X)'
      Do i = 1, nobs
        Write (nout,99999) x(i,1:mx)
      End Do

!     Clean-up the G22 handles
      ifail = 0
      Call g22zaf(hform,ifail)
      Call g22zaf(hddesc,ifail)
      Call g22zaf(hxdesc,ifail)

      Deallocate (dat,x,y,levels,vnames)

99999 Format (100(1X,F4.1))
    End Program g22yafe