Seems the last email didn't deliver the attachement (even though I =20
received it back?)
The attachment is here again but just in case, here's the code as =20
plain text:
#!/usr/bin/env python
#
# ported fnnls.m to fnnls.py
#
# gjok  20050816
#
import sys, math
import numpy, numpy.linalg.linalg as la
def any(X) : return len(numpy.nonzero(X)) !=3D 0
def find(X) : return numpy.nonzero(X)
def norm(X, d) : return max(numpy.sum(abs(X)))
#
# x, w =3D fnnls(XtX, Xty, tol)
#
def fnnls(XtX, Xty, tol =3D 0) :
#FNNLS Nonnegative leastsquares.
#
# Adapted from NNLS of Mathworks, Inc.
# [x,w] =3D nnls(X, y)
#
# x, w =3D fnnls(XtX,Xty) returns the vector X that solves x =3D =
pinv(XtX)=20
*Xty
# in a least squares sense, subject to x >=3D 0.
# Differently stated it solves the problem min y  Xx if
# XtX =3D X'*X and Xty =3D X'*y.
#
# A default tolerance of TOL =3D MAX(SIZE(XtX)) * NORM(XtX,1) * =
EPS
# is used for deciding when elements of x are less than zero.
# This can be overridden with x =3D fnnls(XtX,Xty,TOL).
#
# [x,w] =3D fnnls(XtX,Xty) also returns dual vector w where
# w(i) < 0 where x(i) =3D 0 and w(i) =3D 0 where x(i) > 0.
#
# See also NNLS and FNNLSb
# L. Shure 5887
# Revised, 121588,83189 LS.
# (Partly) Copyright (c) 198494 by The MathWorks, Inc.
# Modified by R. Bro 5796 according to
# Bro R., de Jong S., Journal of Chemometrics, 1997, 11, 393401
# Corresponds to the FNNLSa algorithm in the paper
#=09
# Rasmus bro
# Chemometrics Group, Food Technology
# Dept. Dairy and Food Science
# Royal Vet. & Agricultural
# DK1958 Frederiksberg C
# Denmark
# rb@...
# http://newton.foodsci.kvl.dk/users/rasmus.html
# Reference:
# Lawson and Hanson, "Solving Least Squares Problems", Prentice=20
Hall, 1974.
#
# initialize variables
m =3D XtX.shape[0]
n =3D XtX.shape[1]
if tol =3D=3D 0 :
eps =3D 2.2204e16
tol =3D 10 * eps * norm(XtX,1)*max(m, n);
#end
P =3D numpy.zeros(n, numpy.Int16)
P[:] =3D 1
Z =3D numpy.arange(0,n)
z =3D numpy.zeros(m, numpy.Float)
x =3D numpy.array(P)
ZZ =3D numpy.array(Z)
w =3D Xty  numpy.dot(XtX, x)
# set up iteration criterion
iter =3D 0
itmax =3D 30 * n
# outer loop to put variables into set to hold positive coefficients
while any(Z) and any(w[ZZ] > tol) :
wt =3D w[ZZ].max()
t =3D find(w[ZZ] =3D=3D wt)
t =3D t[1:][0]
t =3D ZZ[t]
P[t] =3D t
Z[t] =3D 1
PP =3D find(P !=3D 1)
ZZ =3D find(Z !=3D 1)
if len(PP) =3D=3D 1 :
XtyPP =3D Xty[PP]
XtXPP =3D XtX[PP, PP]
z[PP] =3D XtyPP / XtXPP
else :
XtyPP =3D numpy.array(Xty[PP])
XtXPP =3D numpy.array(XtX[PP, numpy.array(PP)[:, =20
numpy.NewAxis]])
z[PP] =3D numpy.dot(XtyPP, la.generalized_inverse(XtXPP))
#end
z[ZZ] =3D 0
# inner loop to remove elements from the positive set which no longer =20=
belong
while any(z[PP] <=3D tol) and (iter < itmax) :
iter +=3D 1
iztol =3D find(z <=3D tol)
ip =3D find(P[iztol] !=3D 1)
QQ =3D iztol[ip]
if len(QQ) =3D=3D 1 : alpha =3D x[QQ] / (x[QQ]  z[QQ])
else :
x_xz =3D x[QQ] / (x[QQ]  z[QQ])
alpha =3D x_xz.min()
#end
x +=3D alpha * (z  x)
iabs =3D find(abs(x) < tol)
ip =3D find(P[iabs] !=3D 1)
ij =3D iabs[ip]
Z[ij] =3D numpy.array(ij)
P[ij] =3D 1
PP =3D find(P !=3D 1)
ZZ =3D find(Z !=3D 1)
if len(PP) =3D=3D 1 :
XtyPP =3D Xty[PP]
XtXPP =3D XtX[PP, PP]
z[PP] =3D XtyPP / XtXPP
else :
XtyPP =3D numpy.array(Xty[PP])
XtXPP =3D numpy.array(XtX[PP, numpy.array(PP)[:, =20
numpy.NewAxis]])
z[PP] =3D numpy.dot(XtyPP, =
la.generalized_inverse(XtXPP))
#endif
z[ZZ] =3D 0
#end while
x =3D numpy.array(z)
w =3D Xty  numpy.dot(XtX, x)
#end while
return x, w
#end def
if __name__ =3D=3D '__main__' :
#
# test [x, w] =3D fnnls(Xt.X, Xt.y, tol)
# to solve min y  X.x s.t. x >=3D 0
#
# matlab:lsqnonneg
# X =3D [1, 10, 4, 10; 4, 5, 1, 12; 5, 1, 9, 20];
# y =3D [4; 7; 4]
# x =3D lsqnonneg(X, y) =3D> x =3D [0.9312; 0.3683; 0; 0];
#
X =3D numpy.array([[1, 10, 4, 10], [4, 5, 1, 12], [5, 1, 9, 20]], =20=
numpy.Float)
y =3D numpy.array([4, 7, 4], numpy.Float)
if False :
X =3D numpy.array([[1, 10, 4, 10],
[4, 5, 1, 12],
[5, 1, 9, 20],
[4, 3, 2, 1]], numpy.Float)
y =3D numpy.array([4, 7, 4, 1], numpy.Float)
#end
if False :
X =3D numpy.zeros((20, 20), numpy.Float)
for n in range(20) : X[n,:] =3D numpy.arange(0.0, 400.0, step =20=
=3D 20)
y =3D numpy.arange(0.0, 20.0)
#end
Xt =3D numpy.transpose(numpy.array(X))
x, w =3D fnnls(numpy.dot(Xt, X), numpy.dot(Xt, y))
print 'X =3D ', X
print 'y =3D ', y
print 'x =3D ', x
#end if __name__ =3D=3D '__main__'
Begin forwarded message:
> From: Graeme O'Keefe <gjok@...>
> Date: 25 May 2006 8:46:07 PM
> To: matplotlib user list <matplotlibusers@...>
> Subject: [Matplotlibusers] Fwd: python version of nnls/fnnls
>
>
> Dear matplotlibers,
>
> I seem to use Non Negative Least Squares every couple of months and =20=
> the most recent problem I solved reminded me that I meant to post =20
> this last November.
>
> fnnls.py is a port of fnnls.m which was written by Rasmus Bro =20
> (http://newton.foodsci.kvl.dk/users/rasmus.html, code available at =20
> http://www.ub.es/gesq/mcr/als/fnnls.m, and also available on Matlab =20=
> Central).
>
> I've retained Rasmus' comments and translated the matlab to python/=20
> numpy.
>
=EF=BF=BC
>
> I've used this to linearise problems such as:
> fitting Neutron TOF data with a E0 energy spectrum of gaussians =20=
> (with energyresolutions a function of TOF =3D> E): exp(((E  E0) / =20=
> deltaE)^2)
> fitting Cerebral Blood Flow data with a kspectrum set of =20
> exponentials: exp(kt)
> fitting BGO detector pulses (to detect/correct pileup) with =20
> Tdelayspectrum of exponentials: exp((t  Tdelay) / tau)
>
> Hope others find it as useful as I have.
>
> regards,
>
>
> Graeme
>
>
>
> Graeme O'Keefe, PhD, MACPSEM
> Principal Medical Physicist
> Centre for PET
> Austin Hospital
> Tel: (613)94965767
> Fax: (613) 94585023
>
>
