find_sub.f90

The module LSQ is for unconstrained linear least-squares fitting. It is based upon Applied Statisti

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