! 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