! H05AAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module h05aafe_mod ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Derived Type Definitions .. Type :: calc_subset_data Integer :: n, ldq, ldx Real (Kind=nag_wp) :: tol Real (Kind=nag_wp), Allocatable :: x(:,:), y(:), b(:), cov(:), & h(:), p(:), q(:,:), res(:), & se(:), wk(:), wt(:) Integer, Allocatable :: isx(:) Character (1) :: mean, weight End Type calc_subset_data Contains Subroutine calc_subset_score(m,drop,lz,z,la,a,bscore,cs) ! Calculate the score associated with a particular set of feature subsets. ! M,DROP,LZ,Z,LA,A and BSCORE are all as described in the documentation of ! H05AAF. ! CS - Variable of type CALC_SUBSET_DATA that holds any additional data ! required by this routine ! 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 .. Type (calc_subset_data), Intent (Inout) :: cs Integer, Intent (In) :: drop, la, lz, m ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: bscore(:) Integer, Intent (In) :: a(:), z(:) ! .. Local Scalars .. Real (Kind=nag_wp) :: rss Integer :: i, idf, ifail, inv_drop, ip, & irank Logical :: svd Character (200) :: line ! .. Intrinsic Procedures .. Intrinsic :: abs, allocated, count, max ! .. Executable Statements .. Continue ! Allocate various arrays and read in the data if this is the first time ! this routine has been called If (.Not. allocated(cs%x)) Then ! Read in the number of observations for the data used in the linear ! regression skipping any headings or blank lines Do Read (nin,*,Iostat=ifail) line If (ifail/=0) Exit Read (line,*,Iostat=ifail) cs%n If (ifail==0) Exit End Do ! Read in the control parameters for the linear regression Read (nin,*) cs%tol ! Read in the data cs%ldx = cs%n Allocate (cs%x(cs%ldx,m),cs%y(cs%n)) Read (nin,*)(cs%x(i,1:m),cs%y(i),i=1,cs%n) ! No intercept term and no weights cs%mean = 'Z' cs%weight = 'U' ! Allocate memory required by the regression routine cs%ldq = cs%n Allocate (cs%b(m),cs%cov((m*m+m)/2),cs%h(cs%n),cs%p(m*(m+ & 2)),cs%q(cs%ldq,m+1),cs%res(cs%n),cs%se(m),cs%wk(m*m+5*(m- & 1)),cs%wt(1),cs%isx(m)) End If ! 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) cs%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 cs%isx(z(i)) = inv_drop End Do Do i = 1, max(la,1) If (la>0) Then If (i>1) Then ! Reset the feature altered at the last iteration cs%isx(a(i-1)) = drop End If ! Add or drop the I'th feature in A cs%isx(a(i)) = inv_drop End If ip = count(cs%isx(1:m)==1) ! Fit the regression model ifail = 0 Call g02daf(cs%mean,cs%weight,cs%n,cs%x,cs%ldx,m,cs%isx,ip,cs%y, & cs%wt,rss,idf,cs%b,cs%se,cs%cov,cs%res,cs%h,cs%q,cs%ldq,svd,irank, & cs%p,cs%tol,cs%wk,ifail) ! Return the score (the residual sums of squares) bscore(i) = rss End Do End Subroutine calc_subset_score Subroutine free_subset_score(cs) ! .. Scalar Arguments .. Type (calc_subset_data), Intent (Inout) :: cs ! .. Intrinsic Procedures .. Intrinsic :: allocated ! .. Executable Statements .. If (allocated(cs%x)) Then Deallocate (cs%x) End If If (allocated(cs%y)) Then Deallocate (cs%y) End If If (allocated(cs%b)) Then Deallocate (cs%b) End If If (allocated(cs%cov)) Then Deallocate (cs%cov) End If If (allocated(cs%h)) Then Deallocate (cs%h) End If If (allocated(cs%p)) Then Deallocate (cs%p) End If If (allocated(cs%q)) Then Deallocate (cs%q) End If If (allocated(cs%res)) Then Deallocate (cs%res) End If If (allocated(cs%se)) Then Deallocate (cs%se) End If If (allocated(cs%wk)) Then Deallocate (cs%wk) End If If (allocated(cs%wt)) Then Deallocate (cs%wt) End If If (allocated(cs%isx)) Then Deallocate (cs%isx) End If End Subroutine free_subset_score End Module h05aafe_mod Program h05aafe ! .. Use Statements .. Use nag_library, Only: h05aaf, nag_wp Use h05aafe_mod, Only: calc_subset_data, calc_subset_score, & free_subset_score, nin, nout ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Type (calc_subset_data) :: cs Real (Kind=nag_wp) :: gamma Integer :: cnt, drop, i, ifail, ip, irevcm, & j, la, licomm, lrcomm, lz, m, & mincnt, mincr, mip, nbest ! .. Local Arrays .. Real (Kind=nag_wp) :: acc(2) Real (Kind=nag_wp), Allocatable :: bscore(:), rcomm(:) Integer, Allocatable :: a(:), bz(:,:), ibz(:), icomm(:), & z(:) Logical, Allocatable :: mask(:) ! .. Intrinsic Procedures .. Intrinsic :: max, pack ! .. Executable Statements .. Write (nout,*) 'H05AAF 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))) ! Allocate dummy communication arrays, as we will query the required size licomm = 2 lrcomm = 0 Allocate (icomm(licomm),rcomm(lrcomm)) ! Query the required length for the communication arrays irevcm = 0 ifail = 1 Call h05aaf(irevcm,mincr,m,ip,nbest,drop,lz,z,la,a,bscore,bz,mincnt, & gamma,acc,icomm,licomm,rcomm,lrcomm,ifail) ! Extract the required sizes from the communication array licomm = icomm(1) lrcomm = icomm(2) ! Reallocate communication arrays to the correct size Deallocate (icomm,rcomm) Allocate (icomm(licomm),rcomm(lrcomm)) ! Initialise reverse communication control flag irevcm = 0 ! Call the reverse communication best subset routine in a loop. ! The loop should terminate when IRECVM = 0 on exit cnt = 0 rev_lp: Do ifail = -1 Call h05aaf(irevcm,mincr,m,ip,nbest,drop,lz,z,la,a,bscore,bz,mincnt, & gamma,acc,icomm,licomm,rcomm,lrcomm,ifail) If (irevcm==0) Then If (ifail==0 .Or. ifail==53) Then ! No error, or a warning was issued Exit rev_lp Else ! An error occurred Go To 100 End If End If ! Calculate and return the score for the required models and keep track ! of the number of subsets evaluated cnt = cnt + max(1,la) Call calc_subset_score(m,drop,lz,z,la,a,bscore,cs) End Do rev_lp Call free_subset_score(cs) ! 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==53) Then Write (nout,99997) nbest, & ' subsets of the requested size do not exist, only ', la, & ' are displayed.' End If Write (nout,99996) cnt, ' 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 h05aafe