SUBROUTINE nonnegls(nin, b)
! Non-negative least squares fitting after INCLUD has been used to form
! an orthogonal reduction of the data.
! Author: Alan.Miller @ vic.cmis.csiro.au
! Latest version - 19 October 1998
USE lsq
IMPLICIT NONE
INTEGER, INTENT(OUT) :: nin
REAL (lsq_kind), INTENT(OUT) :: b(:)
! Local variables
REAL (lsq_kind) :: alpha, sxy(ncol), sxx(ncol), temp, wmax, z(ncol)
INTEGER :: i, ier, t
REAL (lsq_kind), PARAMETER :: zero = 0.0_lsq_kind
! The step numbers below are those for Lawson & Hanson's NNLS algorithm
! from page 161 of their book:
! Lawson, C.L. & Hanson, R.J. `Solving least squares problems', Prentice-Hall,
! 1974, republished in SIAM Classics in Applied Mathematics 1995.
! Step 1
nin = 0
b(:ncol) = zero
! Step 2
2 CALL sums_of_sq_and_prod(nin+1, ncol, sxy, sxx)
! Step 4
! We find the max of (sxy)^2 / sxx conditional upon sxy > 0.
! This is better than L & H's (sxy - b'.sxx) which is dependent upon the size
! of the X-variables and so tends to pick X-variables which have large values
! ahead of X-variables with small values.
t = 0
wmax = zero
DO i = nin+1, ncol
IF (sxy(i) > zero) THEN
temp = sxy(i)**2 / sxx(i)
IF (temp > wmax) THEN
wmax = temp
t = i
END IF
END IF
END DO
! Step 3
IF (t == 0) RETURN
! Step 5, add the variable from position t to the model.
nin = nin + 1
IF (t > nin) CALL vmove(t, nin, ier)
! Step 6, calculate the regression coefficients.
6 CALL regcf(z, nin, ier)
! Step 7, if all coefficients > 0, set b = z, then go back to step 2, unless
! all the variables have been included.
IF (ALL(z(:nin) > zero)) THEN
b(:nin) = z(:nin)
IF (nin == ncol) RETURN
GO TO 2
END IF
! Steps 8-11, one or more negative regression coefficients.
! Steps 8-9, calculate shrinkage factor alpha.
t = 0
alpha = 1.0
nin = nin - 1
DO i = 1, nin
IF (z(i) <= zero) THEN
temp = b(i) / (b(i) - z(i))
IF (temp < alpha) THEN
alpha = temp
t = i
END IF
END IF
END DO
! Step 10, shrink regression coefficients
IF (t > 0) THEN
b(1:nin) = b(1:nin) + alpha * (z(1:nin) - b(1:nin))
ELSE
GO TO 2 ! This should never happen
END IF
! Step 11, remove variables with coefficients <= 0.
CALL vmove(t, nin, ier)
nin = nin - 1
GO TO 6
CONTAINS
SUBROUTINE sums_of_sq_and_prod(first, last, sxy, sxx)
! Calculate the sums of squares, and cross-products, of those parts of the
! X and Y variables which are orthogonal to the X-variables already in the
! model (if any) in positions 1 .. first-1. Ignore any variables after last.
INTEGER, INTENT(IN) :: first, last
REAL (lsq_kind), INTENT(OUT) :: sxy(:), sxx(:)
! Local variables
INTEGER :: inc, pos, row, col
REAL (lsq_kind) :: zero = 0.0, diag, dy
! Accumulate sums of squares & products from row FIRST.
! The orthogonal reduction gave:
! X = Q.sqrt(D).R where Q is orthonormal and D is diagonal and such
! that the diagonal elements of R are 1's.
! The sums of squares and products are given by R'DR
! where only rows & columns FIRST .. LAST of R & D are used.
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
RETURN
END SUBROUTINE sums_of_sq_and_prod
END SUBROUTINE nonnegls