!   G22YCF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.

    Module g22ycfe_mod
!     G22YCF 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                           :: print_x, read_line
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6

    Contains
      Subroutine read_line(ierr,v1,v2)
!       Read in a line from NIN, remove any comments and optionally
!       split out the first word

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

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

          If (present(v2)) Then
!           split the first word from the line
            pend = index(v1,' ')
            If (pend/=0) Then
              v2 = adjustl(v1(pend:))
              v1 = adjustl(v1(1:pend))
            Else
              v2 = ''
            End If
          End If
        End If

        Return
      End Subroutine read_line
      Subroutine print_x(intcpt,plab,nobs,mx,x,text)
!       Print the transpose of the first 10 rows of the design matrix

!       .. Parameters ..
        Integer, Parameter             :: max_rows = 10
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: mx, nobs
        Character (*), Intent (In)     :: intcpt, text
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x(:,:)
        Character (*), Intent (In)     :: plab(:)
!       .. Local Scalars ..
        Integer                        :: i, pnobs, si
!       .. Intrinsic Procedures ..
        Intrinsic                      :: min, repeat, trim
!       .. Executable Statements ..
        Continue

!       PLAB holds the labels for the model parameters, so includes a label
!       for the mean effect, if one is present. As the mean effect is not
!       being explicitly included in the design matrix, we may need to skip
!       the first element of PLAB (which will always be the label for the
!       mean effect if one is present)
        If (intcpt=='M') Then
          si = 1
        Else
          si = 0
        End If

!       Printing the first MAX_ROWS rows of the design matrix
        pnobs = min(max_rows,nobs)

!       Display the design matrix
        Write (nout,99998) 'Transpose of First ', pnobs, ' Rows of the ',      &
          text, ' Design Matrix (X)'
        Write (nout,99997) 'Column Name', (i,i=1,pnobs)
        Write (nout,99996) repeat('-',15+pnobs*5)
        Do i = 1, mx
          Write (nout,99999) plab(i+si), x(1:pnobs,i)
        End Do
        Write (nout,*) 'Intercept flag = ', trim(intcpt)

        Return
99999   Format (1X,A15,100(1X,F4.1))
99998   Format (1X,A,I0,A,A,A)
99997   Format (1X,A,3X,100(3X,I2))
99996   Format (1X,A)
      End Subroutine print_x
    End Module g22ycfe_mod

    Program g22ycfe

!     .. Use Statements ..
      Use g22ycfe_mod, Only: nin, nout, print_x, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g22yaf, g22ybf, g22ycf, g22ydf, g22zaf, g22zmf,   &
                             nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hform, hxdesc
      Integer                          :: i, ierr, ifail, ip, lddat, ldx,      &
                                          lisx, lplab, lvinfo, lvnames, mx,    &
                                          nobs, nvar, sddat, sdx
      Character (200)                  :: formula, intcpt, line, tcontrast,    &
                                          tvname
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dat(:,:), x(:,:)
      Real (Kind=nag_wp)               :: tx(0,0)
      Integer, Allocatable             :: isx(:), levels(:), vinfo(:)
      Character (50), Allocatable      :: plab(:), vnames(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: trim
!     .. Executable Statements ..
      Write (nout,*) 'G22YCF 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 formula for the first specification of the model,
!     remove comments and parse it
      Call read_line(ierr,formula)
      ifail = 0
      Call g22yaf(hform,formula,ifail)

!     Read in the contrast to use for all parameters, remove comments and
!     set the contrast optional parameter
      Call read_line(ierr,tcontrast)
      line = 'Contrast=' // trim(tcontrast)
      ifail = 0
      Call g22zmf(hform,line,ifail)

!     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))
      Read (nin,*)(dat(i,1:nvar),i=1,nobs)

!     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/=91) Then
!       redisplay any error messages, other than IFAIL = 91
        ifail = 0
        Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,tx,ldx,sdx,mx,ifail)
      End If

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

!     Generate labels for the columns of X
      lplab = mx + 1
      lvinfo = 0
      lisx = 0
      Allocate (isx(lisx),vinfo(lvinfo),plab(lplab))
      ifail = 0
      Call g22ydf(hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo,vinfo,     &
        ifail)

!     Display the design matrix
      Call print_x(intcpt,plab,nobs,mx,x,'First')

c_lp: Do
!       Read in the name of the variable whose contrasts need to be changed,
!       the value to change them to, remove comments and set the contrast
!       optional argument for the specified variable
        Call read_line(ierr,tvname,tcontrast)
        If (ierr/=0) Then
          Exit c_lp
        End If
        line = 'Contrast:' // trim(tvname) // '=' // trim(tcontrast)
        ifail = 0
        Call g22zmf(hform,line,ifail)
      End Do c_lp

!     Regenerate the design matrix using the new contrasts
!     (the size of X should be the same as previously)
      ifail = 0
      Call g22ycf(hform,hddesc,dat,lddat,sddat,hxdesc,x,ldx,sdx,mx,ifail)

!     Generate labels for the columns of X
      ifail = 0
      Call g22ydf(hform,hxdesc,intcpt,ip,lisx,isx,lplab,plab,lvinfo,vinfo,     &
        ifail)

!     Display the design matrix
      Write (nout,*)
      Call print_x(intcpt,plab,nobs,mx,x,'Second')

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

      Deallocate (dat,x)
      Deallocate (isx,levels,vinfo)
      Deallocate (plab,vnames)

    End Program g22ycfe