MODULE find_subsets
! A module of routines for finding and recording best-fitting subsets of
! regression variables
! Version 1.10, 17 February 2004
! Author: Alan Miller
! formerly of CSIRO Division of Mathematical & Information Sciences
! Phone: (+61) 3 9592-5085
! e-mail: amiller @ bigpond.net.au
! WWW-page: http://users.bigpond.net.au/amiller/
! 17 Feb 2004 Correction to subroutine EFROYM for the case in which all the
! variables are selected. Thanks to David Jones.
! 12 Nov 1999 Made changes to routines exadd1 & add2 to prevent the calculation
! of negative residual sums of squares which could occur in cases
! in which the true RSS is zero. Routine seq2 changed to avoid
! cycling.
! 24 May 2000 Changed lsq_kind to dp (cosmetic change only)
! 4 June 2000 Added routine random_pick which picks a random set of variables.
! 29 Aug 2002 Set value of size in subroutine EFROYM when IER /= 0.
USE lsq
IMPLICIT NONE
INTEGER, SAVE :: max_size, nbest, lopt_dim1
REAL (dp), ALLOCATABLE, SAVE :: bound(:), ress(:,:)
INTEGER, ALLOCATABLE, SAVE :: lopt(:,:)
CONTAINS
SUBROUTINE init_subsets(nvar_max, fit_const, nvar)
INTEGER, INTENT(IN) :: nvar_max
LOGICAL, INTENT(IN) :: fit_const
INTEGER, INTENT(IN), OPTIONAL :: nvar
! Local variables
INTEGER :: i, ier
REAL (dp) :: eps = 1.E-14
LOGICAL :: lindep(ncol)
! The LSQ module has probably already been initialized, but just in case ..
IF (.NOT. initialized) THEN
IF (PRESENT(nvar)) CALL startup(nvar, fit_const)
END IF
IF (fit_const) THEN
max_size = nvar_max + 1
ELSE
max_size = nvar_max
END IF
lopt_dim1 = max_size * (max_size + 1) / 2
IF (ALLOCATED(bound)) DEALLOCATE(bound, ress, lopt)
ALLOCATE (bound(max_size), ress(max_size,nbest), lopt(lopt_dim1,nbest))
bound = HUGE(eps)
ress = HUGE(eps)
lopt = 0
CALL tolset(eps)
CALL sing(lindep, ier)
CALL ss()
DO i = 1, max_size
CALL report(i, rss(i))
END DO
RETURN
END SUBROUTINE init_subsets
SUBROUTINE add1(first, last, ss, smax, jmax, ier)
! Calculate the reduction in residual sum of squares when one variable,
! selected from those in positions FIRST .. LAST, is added in position FIRST,
! given that the variables in positions 1 .. FIRST-1 (if any) are already
! included.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: jmax, ier
REAL (dp), INTENT(OUT) :: ss(:), smax
! Local variables
INTEGER :: j, inc, pos, row, col
REAL (dp) :: zero = 0.0_dp, diag, dy, ssqx, sxx(ncol), sxy(ncol)
! Check call arguments
jmax = 0
smax = zero
ier = 0
IF (first > ncol) ier = 1
IF (last < first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN
! Accumulate sums of squares & products from row FIRST
sxx(first:last) = zero
sxy(first:last) = zero
inc = ncol - last
pos = row_ptr(first)
DO row = first, last
diag = d(row)
dy = diag * rhs(row)
sxx(row) = sxx(row) + diag
sxy(row) = sxy(row) + dy
DO col = row+1, last
sxx(col) = sxx(col) + diag * r(pos)**2
sxy(col) = sxy(col) + dy * r(pos)
pos = pos + 1
END DO
pos = pos + inc
END DO
! Incremental sum of squares for a variable = SXY * SXY / SXX.
! Calculate whenever sqrt(SXX) > TOL for that variable.
DO j = first, last
ssqx = sxx(j)
IF (SQRT(ssqx) > tol(j)) THEN
ss(j) = sxy(j)**2 / sxx(j)
IF (ss(j) > smax) THEN
smax = ss(j)
jmax = j
END IF
ELSE
ss(j) = zero
END IF
END DO
RETURN
END SUBROUTINE add1
SUBROUTINE add2(first, last, smax, j1, j2, ier)
! Calculate the maximum reduction in residual sum of squares when 2
! variables, selected from those in positions FIRST .. LAST, are
! added, given that the variables in positions 1 .. FIRST-1 (if
! any) are already included. J1, J2 are the positions of the two
! best variables. N.B. J2 < J1.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: j1, j2, ier
REAL (dp), INTENT(OUT) :: smax
! Local variables
INTEGER :: start, i1, i2, row, pos1, pos2, inc
REAL (dp) :: zero = 0.0_dp, temp, det, two = 2.0, sxx(ncol), sxy(ncol), sx1x2
! Check call arguments
smax = zero
j1 = 0
j2 = 0
ier = 0
IF (first > ncol) ier = 1
IF (last <= first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN
start = row_ptr(first)
! Cycle through all pairs of variables from those between FIRST & LAST.
DO i1 = first, last
sxx(i1) = d(i1)
sxy(i1) = d(i1) * rhs(i1)
pos1 = start + i1 - first - 1
DO row = first, i1-1
temp = d(row) * r(pos1)
sxx(i1) = sxx(i1) + temp*r(pos1)
sxy(i1) = sxy(i1) + temp*rhs(row)
pos1 = pos1 + ncol - row - 1
END DO
DO i2 = first, i1-1
pos1 = start + i1 - first - 1
pos2 = start + i2 - first - 1
sx1x2 = zero
DO row = first, i2-1
sx1x2 = sx1x2 + d(row)*r(pos1)*r(pos2)
inc = ncol - row - 1
pos1 = pos1 + inc
pos2 = pos2 + inc
END DO
sx1x2 = sx1x2 + d(i2)*r(pos1)
! Calculate reduction in RSS for pair I1, I2.
! The sum of squares & cross-products are in:
! ( SXX(I1) SX1X2 ) ( SXY(I1) )
! ( SX1X2 SXX(I2) ) ( SXY(I2) )
det = MAX( (sxx(i1) * sxx(i2) - sx1x2**2), zero)
temp = SQRT(det)
IF (temp < tol(i1)*SQRT(sxx(i2)) .OR. &
temp < tol(i2)*SQRT(sxx(i1))) CYCLE
temp = ((sxx(i2)*sxy(i1) - two*sx1x2*sxy(i2))*sxy(i1) + sxx(i1)*sxy(i2)**2) &
/ det
IF (temp > smax) THEN
smax = temp
j1 = i1
j2 = i2
END IF
END DO ! i2 = first, i1-1
END DO ! i1 = first, last
RETURN
END SUBROUTINE add2
SUBROUTINE bakwrd(first, last, ier)
! Backward elimination from variables in positions FIRST .. LAST.
! If FIRST > 1, variables in positions prior to this are forced in.
! If LAST < ncol, variables in positions after this are forced out.
! On exit, the array VORDER contains the numbers of the variables
! in the order in which they were deleted.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: ier
! Local variables
INTEGER :: pos, jmin, i
REAL (dp) :: ss(last), smin
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN
! For POS = LAST, ..., FIRST+1 call DROP1 to find best variable to
! find which variable to drop next.
DO pos = last, first+1, -1
CALL drop1(first, pos, ss, smin, jmin, ier)
CALL exdrop1(first, pos, ss, smin, jmin)
IF (jmin > 0 .AND. jmin < pos) THEN
CALL vmove(jmin, pos, ier)
IF (nbest > 0) THEN
DO i = jmin, pos-1
CALL report(i, rss(i))
END DO
END IF
END IF
END DO
RETURN
END SUBROUTINE bakwrd
SUBROUTINE drop1(first, last, ss, smin, jmin, ier)
! Calculate the increase in the residual sum of squares when the variable in
! position J is dropped from the model (i.e. moved to position LAST),
! for J = FIRST, ..., LAST-1.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: jmin, ier
REAL (dp), INTENT(OUT) :: ss(:), smin
! Local variables
INTEGER :: j, pos1, inc, pos, row, col, i
REAL (dp) :: large = HUGE(1.0_dp), zero = 0.0_dp, d1, rhs1, d2, x, wk(last), &
vsmall = TINY(1.0_dp)
! Check call arguments
jmin = 0
smin = large
ier = 0
IF (first > ncol) ier = 1
IF (last < first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN
! POS1 = position of first element of row FIRST in r.
pos1 = row_ptr(first)
inc = ncol - last
! Start of outer cycle for the variable to be dropped.
DO j = first, last
d1 = d(j)
IF (SQRT(d1) < tol(j)) THEN
ss(j) = zero
smin = zero
jmin = j
GO TO 50
END IF
rhs1 = rhs(j)
IF (j == last) GO TO 40
! Copy row J of R into WK.
pos = pos1
DO i = j+1, last
wk(i) = r(pos)
pos = pos + 1
END DO
pos = pos + inc
! Lower the variable past each row.
DO row = j+1, last
x = wk(row)
d2 = d(row)
IF (ABS(x) * SQRT(d1) < tol(row) .OR. d2 < vsmall) THEN
pos = pos + ncol - row
CYCLE
END IF
d1 = d1 * d2 / (d2 + d1 * x**2)
DO col = row+1, last
wk(col) = wk(col) - x * r(pos)
pos = pos + 1
END DO
rhs1 = rhs1 - x * rhs(row)
pos = pos + inc
END DO
40 ss(j) = rhs1 * d1 * rhs1
IF (ss(j) < smin) THEN
jmin = j
smin = ss(j)
END IF
! Update position of first element in row of r.
50 IF (j < last) pos1 = pos1 + ncol - j
END DO
RETURN
END SUBROUTINE drop1
SUBROUTINE efroym(first, last, fin, fout, size, ier, lout)
! Efroymson's stepwise regression from variables in positions FIRST,
! ..., LAST. If FIRST > 1, variables in positions prior to this are
! forced in. If LAST < ncol, variables in positions after this are
! forced out.
! A report is written to unit LOUT if LOUT >= 0.
INTEGER, INTENT(IN) :: first, last, lout
INTEGER, INTENT(OUT) :: size, ier
REAL (dp), INTENT(IN) :: fin, fout
! Local variables
INTEGER :: jmax, jmin, i
REAL (dp) :: one = 1.0, eps, zero = 0.0, ss(last), smax, base, var, f, smin
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (fin < fout .OR. fin <= zero) ier = ier + 256
IF (nobs <= ncol) ier = ier + 512
IF (ier /= 0) THEN
size = 0
RETURN
END IF
! EPS approximates the smallest quantity such that the calculated value of
! (1 + EPS) is > 1. It is used to test for a perfect fit (RSS = 0).
eps = EPSILON(one)
! SIZE = number of variables in the current subset
size = first - 1
! Find the best variable to add next
20 CALL add1(size+1, last, ss, smax, jmax, ier)
IF (nbest > 0) CALL exadd1(size+1, smax, jmax, ss, last)
! Calculate 'F-to-enter' value
IF (size > 0) THEN
base = rss(size)
ELSE
base = rss(1) + ss(1)
END IF
var = (base - smax) / (nobs - size - 1)
IF (var < eps*base) THEN
ier = -1
f = zero
ELSE
f = smax / var
END IF
IF (lout >= 0) WRITE(lout, 900) vorder(jmax), f
900 FORMAT(' Best variable to add: ', i4, ' F-to-enter = ', f10.2)
! Exit if F < FIN or IER < 0 (perfect fit)
IF (f < fin .OR. ier < 0) RETURN
! Add the variable to the subset (in position FIRST).
IF (lout >= 0) WRITE(lout, '(50x, "Variable added")')
size = size + 1
IF (jmax > first) CALL vmove(jmax, first, ier)
DO i = first, MIN(jmax-1, max_size)
CALL report(i, rss(i))
END DO
! See whether a variable entered earlier can be deleted now.
30 IF (size <= first) GO TO 20
CALL drop1(first+1, size, ss, smin, jmin, ier)
CALL exdrop1(first+1, size, ss, smin, jmin)
var = rss(size) / (nobs - size)
f = smin / var
IF (lout >= 0) WRITE(lout, 910) vorder(jmin), f
910 FORMAT(' Best variable to drop: ', i4, ' F-to-drop = ', f10.2)
IF (f < fout) THEN
IF (lout >= 0) WRITE(lout, '(50x, "Variable dropped")')
CALL vmove(jmin, size, ier)
IF (nbest > 0) THEN
DO i = jmin, size-1
CALL report(i, rss(i))
END DO
END IF
size = size - 1
GO TO 30
END IF
IF (size >= last) RETURN
GO TO 20
END SUBROUTINE efroym
SUBROUTINE exadd1(ivar, smax, jmax, ss, last)
! Update the NBEST subsets of IVAR variables found from a call
! to subroutine ADD1.
INTEGER, INTENT(IN) :: ivar, jmax, last
REAL (dp), INTENT(IN) :: smax, ss(:)
! Local variables
REAL (dp) :: zero = 0.0_dp, ssbase, sm, temp, wk(last)
INTEGER :: i, j, ltemp, jm
IF (jmax == 0) RETURN
IF (ivar <= 0) RETURN
IF (ivar > max_size) RETURN
ltemp = vorder(ivar)
jm = jmax
sm = smax
IF (ivar > 1) ssbase = rss(ivar-1)
IF (ivar == 1) ssbase = rss(ivar) + ss(1)
wk(ivar:last) = ss(ivar:last)
DO i = 1, nbest
temp = MAX(ssbase - sm, zero)
IF (temp >= bound(ivar)) EXIT
vorder(ivar) = vorder(jm)
IF (jm == ivar) vorder(ivar) = ltemp
CALL report(ivar, temp)
IF (i >= nbest) EXIT
wk(jm) = zero
sm = zero
jm = 0
DO j = ivar, last
IF (wk(j) <= sm) CYCLE
jm = j
sm = wk(j)
END DO
IF (jm == 0) EXIT
END DO
! Restore VORDER(IVAR)
vorder(ivar) = ltemp
RETURN
END SUBROUTINE exadd1
SUBROUTINE exdrop1(first, last, ss, smin, jmin)
! Record any new subsets of (LAST-1) variables found from a call to DROP1
INTEGER, INTENT(IN) :: first, last, jmin
REAL (dp), INTENT(IN) :: ss(:), smin
! Local variables
INTEGER :: list(1:last), i
REAL (dp) :: rss_last, ssq
IF (jmin == 0 .OR. last < 1 .OR. last-1 > max_size) RETURN
rss_last = rss(last)
IF (rss_last + smin > bound(last-1)) RETURN
list = vorder(1:last)
DO i = first, last-1
vorder(i:last-1) = list(i+1:last)
ssq = rss_last + ss(i)
CALL report(last-1, ssq)
vorder(i) = list(i)
END DO
RETURN
END SUBROUTINE exdrop1
SUBROUTINE forwrd(first, last, ier)
! Forward selection from variables in positions FIRST .. LAST.
! If FIRST > 1, variables in positions prior to this are forced in.
! If LAST < ncol, variables in positions after this are forced out.
! On exit, the array VORDER contains the numbers of the variables
! in the order in which they were added.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: ier
! Local variables
INTEGER :: pos, jmax
REAL (dp) :: ss(last), smax
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN
! For POS = FIRST .. max_size, call ADD1 to find best variable to put
! into position POS.
DO pos = first, max_size
CALL add1(pos, last, ss, smax, jmax, ier)
IF (nbest > 0) CALL exadd1(pos, smax, jmax, ss, last)
! Move the best variable to position POS.
IF (jmax > pos) CALL vmove(jmax, pos, ier)
END DO
RETURN
END SUBROUTINE forwrd
SUBROUTINE report(nv, ssq)
! Update record of the best NBEST subsets of NV variables, if
! necessary, using SSQ.
INTEGER, INTENT(IN) :: nv
REAL (dp), INTENT(IN) :: ssq
! Local variables
INTEGER :: rank, pos1, j, list(nv)
REAL (dp) :: under1 = 0.99999_dp, above1 = 1.00001_dp
! If residual sum of squares (SSQ) for the new subset > the
! appropriate bound, return.
IF(nv > max_size) RETURN
IF(ssq >= bound(nv)) RETURN
pos1 = (nv*(nv-1))/2 + 1
! Find rank of the new subset
DO rank = 1, nbest
IF(ssq < ress(nv,rank)*above1) THEN
list = vorder(1:nv)
CALL shell(list, nv)
! Check list of variables if ssq is almost equal to ress(nv,rank) -
! to avoid including the same subset twice.
IF (ssq > ress(nv,rank)*under1) THEN
IF (same_vars(list, lopt(pos1:,rank), nv)) RETURN
END IF
! Record the new subset, and move the others down one place.
DO j = nbest-1, rank, -1
ress(nv,j+1) = ress(nv,j)
lopt(pos1:pos1+nv-1, j+1) = lopt(pos1:pos1+nv-1, j)
END DO
ress(nv,rank) = ssq
lopt(pos1:pos1+nv-1, rank) = list(1:nv)
bound(nv) = ress(nv,nbest)
RETURN
END IF
END DO
RETURN
END SUBROUTINE report
SUBROUTINE shell(l, n)
! Perform a SHELL-sort on integer array L, sorting into increasing order.
! Latest revision - 5 July 1995
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN OUT) :: l(:)
! Local variables
INTEGER :: start, finish, temp, new, i1, i2, incr, it
incr = n
DO
incr = incr/3
IF (incr == 2*(incr/2)) incr = incr + 1
DO start = 1, incr
finish = n
! TEMP contains the element being compared; IT holds its current
! location. It is compared with the elements in locations
! IT+INCR, IT+2.INCR, ... until a larger element is found. All
! smaller elements move INCR locations towards the start. After
! each time through the sequence, the FINISH is decreased by INCR
! until FINISH <= INCR.
20 i1 = start
temp = l(i1)
it = i1
! I2 = location of element new to be compared with TEMP.
! Test I2 <= FINISH.
DO
i2 = i1 + incr
IF (i2 > finish) THEN
IF (i1 > it) l(i1) = temp
finish = finish - incr
EXIT
END IF
new = l(i2)
! If TEMP > NEW, move NEW to lower-numbered position.
IF (temp > new) THEN
l(i1) = new
i1 = i2
CYCLE
END IF
! TEMP <= NEW so do not swap.
! Use NEW as the next TEMP.
IF (i1 > it) l(i1) = temp
i1 = i2
temp = new
it = i1
! Repeat until FINISH <= INCR.
END DO
IF (finish > incr) GO TO 20
END DO
! Repeat until INCR = 1.
IF (incr <= 1) RETURN
END DO
RETURN
END SUBROUTINE shell
FUNCTION same_vars(list1, list2, n) RESULT(same)
LOGICAL :: same
INTEGER, INTENT(IN) :: n, list1(:), list2(:)
same = ALL(list1(1:n) == list2(1:n))
RETURN
END FUNCTION same_vars
SUBROUTINE seq2(first, last, ier)
! Sequential replacement algorithm applied to the variables in positions
! FIRST, ..., LAST. 2 variables at a time are added or replaced.
! If FIRST > 1, variables in positions prior to this are forced in.
! If LAST < NP, variables in positions after this are left out.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: ier
! Local variables
INTEGER :: nv, nsize
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0 .OR. nbest <= 0) RETURN
nv = MIN(max_size, last-1)
! Outer loop; SIZE = current size of subset being considered.
DO nsize = first+1, nv
CALL replace2(first, last, nsize)
END DO
RETURN
END SUBROUTINE seq2
SUBROUTINE replace2(first, last, nsize)
! Replace 2 variables at a time from those in positions first, ..., nsize
! with 2 from positions nsize, .., last - if they reduce the RSS.
INTEGER, INTENT(IN) :: first, last, nsize
! Local variables
INTEGER :: ier, j1, j2, pos1, pos2, best(2), i, iwk(last)
REAL (dp) :: smax, rssnew, rssmin, save_rss
REAL (dp), PARAMETER :: zero = 0.0_dp
10 best(1) = 0
best(2) = 0
rssmin = rss(nsize)
! Two loops to place all pairs of variables in positions nsize-1 and nsize.
! POS1 = destination for variable from position nsize.
! POS2 = destination for variable from position nsize-1.
DO pos1 = first, nsize
DO pos2 = pos1, nsize-1
CALL add2(nsize-1, last, smax, j1, j2, ier)
IF (j1+j2 > nsize + nsize - 1) THEN
rssnew = MAX(rss(nsize-2) - smax, zero)
IF (rssnew < rssmin) THEN
best(1) = vorder(j1)
best(2) = vorder(j2)
iwk(1:nsize-2) = vorder(1:nsize-2)
rssmin = rssnew
END IF
END IF
CALL vmove(nsize-1, pos2, ier)
END DO
CALL vmove(nsize, pos1, ier)
DO i = pos1, nsize
CALL report(i, rss(i))
END DO
END DO
! If any replacement reduces the RSS, make the best one.
IF (best(1) + best(2) > 0) THEN
iwk(nsize-1) = best(2)
iwk(nsize) = best(1)
save_rss = rss(nsize)
CALL reordr(iwk, nsize, 1, ier)
DO i = first, nsize
CALL report(i, rss(i))
END DO
! The calculated value of rssmin above is only a rough approximation to
! the real residual sum of squares, thiugh usually good enough.
! The new value of rss(nsize) is more accurate. It is used below
! to avoid cycling when several subsets give the same RSS.
IF (rss(nsize) < save_rss) GO TO 10
END IF
RETURN
END SUBROUTINE replace2
SUBROUTINE seqrep(first, last, ier)
! Sequential replacement algorithm applied to the variables in
! positions FIRST, ..., LAST.
! If FIRST > 1, variables in positions prior to this are forced in.
! If LAST < ncol, variables in positions after this are forced out.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: ier
! Local variables
INTEGER :: nv, size, start, best, from, i, jmax, count, j
REAL (dp) :: zero = 0.0_dp, ssred, ss(last), smax
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0 .OR. nbest <= 0) RETURN
nv = MIN(max_size, last-1)
! Outer loop; SIZE = current size of subset being considered.
DO size = first, nv
count = 0
start = first
10 ssred = zero
best = 0
from = 0
! Find the best variable from those in positions SIZE+1, ..., LAST
! to replace the one in position SIZE. Then rotate variables in
! positions START, ..., SIZE.
DO i = start, size
CALL add1(size, last, ss, smax, jmax, ier)
IF (jmax > size) THEN
CALL exadd1(size, smax, jmax, ss, last)
IF (smax > ssred) THEN
ssred = smax
best = jmax
IF (i < size) THEN
from = size + start - i - 1
ELSE
from = size
END IF
END IF
END IF
IF (i < size) CALL vmove(size, start, ier)
DO j = start, size-1
CALL report(j, rss(j))
END DO
END DO ! i = start, size
! If any replacement reduces the RSS, make the best one.
! Move variable from position FROM to SIZE.
! Move variable from position BEST to FIRST.
IF (best > size) THEN
IF (from < size) CALL vmove(from, size, ier)
CALL vmove(best, first, ier)
DO j = first, best-1
CALL report(j, rss(j))
END DO
count = 0
start = first + 1
ELSE
count = count + 1
END IF
! Repeat until COUNT = SIZE - START + 1
IF (count <= size - start) GO TO 10
END DO
RETURN
END SUBROUTINE seqrep
SUBROUTINE xhaust(first, last, ier)
! Exhaustive search algorithm, using leaps and bounds, applied to
! the variables in positions FIRST, ..., LAST.
! If FIRST > 1, variables in positions prior to this are forced in.
! If LAST < ncol, variables in positions after this are forced out.
INTEGER, INTENT(IN) :: first, last
INTEGER, INTENT(OUT) :: ier
! Local variables
INTEGER :: row, i, jmax, ipt, newpos, iwk(max_size)
REAL (dp) :: ss(last), smax, temp
! Check call arguments
ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0 .OR. nbest <= 0) RETURN
! Record subsets contained in the initial ordering, including check
! for variables which are linearly related to earlier variables.
! This should be redundant if the user has first called SING and
! init_subsets.
DO row = first, max_size
IF (d(row) <= tol(row)) THEN
ier = -999
RETURN
END IF
CALL report(row, rss(row))
END DO
! IWK(I) contains the upper limit for the I-th simulated DO-loop for
! I = FIRST, ..., max_size-1.
! IPT points to the current DO loop.
iwk(first:max_size) = last
! Innermost loop.
! Find best possible variable for position max_size from those in
! positions max_size, .., IWK(max_size).
30 CALL add1(max_size, iwk(max_size), ss, smax, jmax, ier)
CALL exadd1(max_size, smax, jmax, ss, iwk(max_size))
! Move to next lower numbered loop which has not been exhausted.
ipt = max_size - 1
40 IF (ipt >= iwk(ipt)) THEN
ipt = ipt - 1
IF (ipt >= first) GO TO 40
RETURN
END IF
! Lower variable from position IPT to position IWK(IPT).
! Record any good new subsets found by the move.
newpos = iwk(ipt)
CALL vmove(ipt, newpos, ier)
DO i = ipt, MIN(max_size, newpos-1)
CALL report(i, rss(i))
END DO
! Reset all ends of loops for I >= IPT.
iwk(ipt:max_size) = newpos - 1
! If residual sum of squares for all variables above position NEWPOS
! is greater than BOUND(I), no better subsets of size I can be found
! inside the current loop.
temp = rss(newpos-1)
DO i = ipt, max_size
IF (temp > bound(i)) GO TO 80
END DO
IF (iwk(max_size) > max_size) GO TO 30
ipt = max_size - 1
GO TO 40
80 ipt = i - 1
IF (ipt < first) RETURN
GO TO 40
END SUBROUTINE xhaust
SUBROUTINE random_pick(first, last, npick)
! Pick npick variables at random from those in positions first, ..., last
! and move them to occupy positions starting from first.
INTEGER, INTENT(IN) :: first, last, npick
! Local variables
INTEGER :: first2, i, ilist(1:last), j, k, navail
REAL :: r
navail = last + 1 - first
IF (npick >= navail .OR. npick <= 0) RETURN
DO i = first, last
ilist(i) = vorder(i)
END DO
first2 = first
DO i = 1, npick
CALL RANDOM_NUMBER(r)
k = first2 + r*navail
IF (k > first2) THEN
j = ilist(first2)
ilist(first2) = ilist(k)
ilist(k) = j
END IF
first2 = first2 + 1
navail = navail - 1
END DO
CALL reordr(ilist(first:), npick, first, i)
RETURN
END SUBROUTINE random_pick
END MODULE find_subsets