Content-Type: multipart/mixed; boundary=Apple-Mail-27-613260171 --Apple-Mail-27-613260171 Content-Transfer-Encoding: quoted-printable Content-Type: text/html; charset=ISO-8859-1 Seems the last email didn't = deliver the attachement (even though I received it back?)

The attachment is here = again but just in case, here's the code as 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)=A0 =A0=A0 : = return len(numpy.nonzero(X)) !=3D 0
def find(X)=A0 =A0 : 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 Non-negative = least-squares.
#
# Adapted from NNLS of Mathworks, = Inc.
#=A0= =A0 =A0 =A0 =A0 [x,w] =3D nnls(X, y)
#
# x, w =3D fnnls(XtX,Xty) = returns the vector X that solves x =3D = pinv(XtX)*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.
#

# L. Shure = 5-8-87
# Revised, 12-15-88,8-31-89 LS.
# (Partly) Copyright (c) = 1984-94 by The MathWorks, Inc.

# Modified by R. Bro 5-7-96 = according to
#=A0 =A0 =A0=A0 Bro R., de Jong S., Journal of = Chemometrics, 1997, 11, 393-401
# Corresponds to the FNNLSa = algorithm in the paper
#
# Rasmus = bro
# Chemometrics Group, Food = Technology
# Dept. Dairy and Food Science
# Royal Vet. & = Agricultural
# DK-1958 Frederiksberg C
# Denmark
# rb@kvl.dk
#=A0 = Reference:
#=A0 Lawson and Hanson, "Solving Least Squares = Problems", Prentice-Hall, 1974.
#

# initialize = variables
=A0 =A0 m =3D XtX.shape[0]
=A0 =A0 n =3D = XtX.shape[1]

=A0 =A0 if tol =3D=3D 0 = :
=A0 = =A0 =A0 =A0 eps =3D 2.2204e-16
=A0 =A0 =A0 =A0 tol =3D 10 * = eps * norm(XtX,1)*max(m, n);
=A0 =A0 = #end

=A0 =A0 P =3D numpy.zeros(n, = numpy.Int16)
=A0 =A0 P[:] =3D -1
=A0 =A0 Z =3D = numpy.arange(0,n)

=A0 =A0 z =3D numpy.zeros(m, = numpy.Float)
=A0 =A0 x =3D numpy.array(P)
=A0 =A0 ZZ =3D = numpy.array(Z)

=A0 =A0 w =3D Xty - = numpy.dot(XtX, x)

=A0 # set up iteration = criterion
=A0 =A0 iter =3D 0
=A0 =A0 itmax =3D 30 * = n

=A0 # outer loop to put = variables into set to hold positive coefficients
=A0 =A0 while any(Z) and = any(w[ZZ] > tol) :
=A0 =A0 =A0 =A0 wt =3D = w[ZZ].max()
=A0 =A0 =A0 =A0 t =3D find(w[ZZ] =3D=3D = wt)
=A0= =A0 =A0 =A0 t =3D t[-1:][0]
=A0 =A0 =A0 =A0 t =3D = ZZ[t]
=A0= =A0 =A0 =A0 P[t] =3D t
=A0 =A0 =A0 =A0 Z[t] =3D -1
=A0 =A0 =A0 =A0 PP =3D = find(P !=3D -1)

=A0 =A0 =A0 =A0 ZZ =3D = find(Z !=3D -1)
=A0 =A0 =A0 =A0 if len(PP) =3D=3D 1 = :
=A0 = =A0 =A0 =A0 =A0 =A0 XtyPP =3D Xty[PP]
=A0 =A0 =A0 =A0 =A0 =A0 = XtXPP =3D XtX[PP, PP]
=A0 =A0 =A0 =A0 =A0 =A0 z[PP] =3D XtyPP / = XtXPP
=A0= =A0 =A0 =A0 else :
=A0 =A0 =A0 =A0 =A0 =A0 XtyPP =3D = numpy.array(Xty[PP])
=A0 =A0 =A0 =A0 =A0 =A0 XtXPP =3D numpy.array(XtX[PP, = numpy.array(PP)[:, numpy.NewAxis]])
=A0 =A0 =A0 =A0 =A0 =A0 = z[PP] =3D numpy.dot(XtyPP, = la.generalized_inverse(XtXPP))
=A0 =A0 =A0 =A0 = #end
=A0= =A0 =A0 =A0 z[ZZ] =3D 0

# inner loop to remove = elements from the positive set which no longer = belong
=A0 =A0 =A0 =A0 while any(z[PP] <=3D tol) and (iter = < itmax) :
=A0 =A0 =A0 =A0 =A0 =A0 iter +=3D = 1
=A0 = =A0 =A0 =A0 =A0 =A0 iztol =3D find(z <=3D tol)
=A0 =A0 =A0 =A0 =A0 =A0 ip =3D= find(P[iztol] !=3D -1)
=A0 =A0 =A0 =A0 =A0 =A0 QQ =3D = iztol[ip]

=A0 =A0 =A0 =A0 =A0 =A0 if = len(QQ) =3D=3D 1 : alpha =3D x[QQ] / (x[QQ] - = z[QQ])
=A0 =A0 =A0 =A0 =A0 =A0 else :
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 x_xz =3D x[QQ] / (x[QQ] - z[QQ])
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 alpha =3D x_xz.min()
=A0 =A0 =A0 =A0 =A0 =A0 #end

=A0 =A0 =A0 =A0 =A0 =A0 x +=3D= alpha * (z - x)
=A0 =A0 =A0 =A0 =A0 =A0 iabs =3D find(abs(x) < = tol)
=A0= =A0 =A0 =A0 =A0 =A0 ip =3D find(P[iabs] !=3D -1)
=A0 =A0 =A0 =A0 =A0 =A0 ij =3D= iabs[ip]

=A0 =A0 =A0 =A0 =A0 =A0 = Z[ij] =3D numpy.array(ij)
=A0 =A0 =A0 =A0 =A0 =A0 = P[ij] =3D -1
=A0 =A0 =A0 =A0 =A0 =A0 PP =3D find(P !=3D = -1)
=A0= =A0 =A0 =A0 =A0 =A0 ZZ =3D find(Z !=3D -1)

=A0 =A0 =A0 =A0 =A0 =A0 if = len(PP) =3D=3D 1 :
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 XtyPP =3D = Xty[PP]
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 XtXPP =3D XtX[PP, = PP]
=A0= =A0 =A0 =A0 =A0 =A0 =A0 =A0 z[PP] =3D XtyPP / = XtXPP
=A0= =A0 =A0 =A0 =A0 =A0 else :
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 XtyPP =3D numpy.array(Xty[PP])
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 XtXPP =3D numpy.array(XtX[PP, numpy.array(PP)[:, = numpy.NewAxis]])
=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 z[PP] =3D = numpy.dot(XtyPP, la.generalized_inverse(XtXPP))
=A0 =A0 =A0 =A0 =A0 =A0 = #endif
=A0 =A0 =A0 =A0 =A0 =A0 z[ZZ] =3D = 0
=A0 = =A0 =A0 =A0 #end while
=A0 =A0 =A0 =A0 x =3D = numpy.array(z)
=A0 =A0 =A0 =A0 w =3D Xty - numpy.dot(XtX, = x)
=A0 = =A0 #end while

=A0 =A0 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
#=A0=A0 X =3D [1, 10, 4, 10; 4, 5, 1, 12; 5, 1, 9, = 20];
#=A0=A0 y =3D [4; 7; 4]
#=A0=A0 x =3D lsqnonneg(X, = y) =3D> x =3D [0.9312; 0.3683; 0; 0];
#

=A0 =A0 X =3D = numpy.array([[1, 10, 4, 10], [4, 5, 1, 12], [5, 1, 9, 20]], = numpy.Float)
=A0 =A0 y =3D numpy.array([4, 7, 4], = numpy.Float)

=A0 =A0 if False = :
=A0 = =A0 =A0 =A0 X =3D numpy.array([[1, 10, 4, 10],
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 =A0 =A0 =A0 =A0=A0 [4, 5, 1, 12],
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 =A0 =A0 =A0 =A0=A0 [5, 1, 9, 20],
=A0 =A0 =A0 =A0 =A0 =A0 =A0 = =A0 =A0 =A0 =A0 =A0=A0 [4, 3, 2, 1]], = numpy.Float)
=A0 =A0 =A0 =A0 y =3D numpy.array([4, 7, 4, 1], = numpy.Float)
=A0 =A0 #end

=A0 =A0 if False = :
=A0 = =A0 =A0 =A0 X =3D numpy.zeros((20, 20), = numpy.Float)
=A0 =A0 =A0 =A0 for n in range(20) : X[n,:] =3D = numpy.arange(0.0, 400.0, step =3D 20)
=A0 =A0 =A0 =A0 y =3D = numpy.arange(0.0, 20.0)
=A0 =A0 #end
=A0 =A0 Xt =3D = numpy.transpose(numpy.array(X))
=A0 =A0 x, w =3D = fnnls(numpy.dot(Xt, X), numpy.dot(Xt, y))

=A0 =A0 print 'X =3D ', = X
=A0 = =A0 print 'y =3D ', y
=A0 =A0 print 'x =3D ', x
#end if __name__ =3D=3D = '__main__'

Begin forwarded = message:

From: Graeme O'Keefe <gjok@netspace.net.au>
Date: 25 May 2006 8:46:07 PM
To: matplotlib user list <matplotlib-users@li= sts.sourceforge.net>
Subject: [Matplotlib-users] Fwd: = python version of nnls/fnnls

Dear matplotlibers,

I seem to use Non Negative = Least Squares every couple of months and the most recent problem I = solved reminded me that I meant to post this last = November.

fnnls.py is a port of fnnls.m which was = written by Rasmus Bro (http://newton.food= sci.kvl.dk/users/rasmus.html, code available at=A0http://www.ub.es/gesq/mcr/a= ls/fnnls.m, and also available on Matlab = Central).

I've = retained Rasmus' comments and translated the matlab to = python/numpy.

= --Apple-Mail-27-613260171 Content-Transfer-Encoding: 7bit Content-Type: text/x-python-script; x-unix-mode=0755; x-mac-creator=454D4178; name="fnnls.py" Content-Disposition: attachment; filename=fnnls.py #!/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)) != 0 def find(X) : return numpy.nonzero(X) def norm(X, d) : return max(numpy.sum(abs(X))) # # x, w = fnnls(XtX, Xty, tol) # def fnnls(XtX, Xty, tol = 0) : #FNNLS Non-negative least-squares. # # Adapted from NNLS of Mathworks, Inc. # [x,w] = nnls(X, y) # # x, w = fnnls(XtX,Xty) returns the vector X that solves x = pinv(XtX)*Xty # in a least squares sense, subject to x >= 0. # Differently stated it solves the problem min ||y - Xx|| if # XtX = X'*X and Xty = X'*y. # # A default tolerance of TOL = 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 = fnnls(XtX,Xty,TOL). # # [x,w] = fnnls(XtX,Xty) also returns dual vector w where # w(i) < 0 where x(i) = 0 and w(i) = 0 where x(i) > 0. # # See also NNLS and FNNLSb # L. Shure 5-8-87 # Revised, 12-15-88,8-31-89 LS. # (Partly) Copyright (c) 1984-94 by The MathWorks, Inc. # Modified by R. Bro 5-7-96 according to # Bro R., de Jong S., Journal of Chemometrics, 1997, 11, 393-401 # Corresponds to the FNNLSa algorithm in the paper # # Rasmus bro # Chemometrics Group, Food Technology # Dept. Dairy and Food Science # Royal Vet. & Agricultural # DK-1958 Frederiksberg C # Denmark # rb@kvl.dk # http://newton.foodsci.kvl.dk/users/rasmus.html # Reference: # Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974. # # initialize variables m = XtX.shape[0] n = XtX.shape[1] if tol == 0 : eps = 2.2204e-16 tol = 10 * eps * norm(XtX,1)*max(m, n); #end P = numpy.zeros(n, numpy.Int16) P[:] = -1 Z = numpy.arange(0,n) z = numpy.zeros(m, numpy.Float) x = numpy.array(P) ZZ = numpy.array(Z) w = Xty - numpy.dot(XtX, x) # set up iteration criterion iter = 0 itmax = 30 * n # outer loop to put variables into set to hold positive coefficients while any(Z) and any(w[ZZ] > tol) : wt = w[ZZ].max() t = find(w[ZZ] == wt) t = t[-1:][0] t = ZZ[t] P[t] = t Z[t] = -1 PP = find(P != -1) ZZ = find(Z != -1) if len(PP) == 1 : XtyPP = Xty[PP] XtXPP = XtX[PP, PP] z[PP] = XtyPP / XtXPP else : XtyPP = numpy.array(Xty[PP]) XtXPP = numpy.array(XtX[PP, numpy.array(PP)[:, numpy.NewAxis]]) z[PP] = numpy.dot(XtyPP, la.generalized_inverse(XtXPP)) #end z[ZZ] = 0 # inner loop to remove elements from the positive set which no longer belong while any(z[PP] <= tol) and (iter < itmax) : iter += 1 iztol = find(z <= tol) ip = find(P[iztol] != -1) QQ = iztol[ip] if len(QQ) == 1 : alpha = x[QQ] / (x[QQ] - z[QQ]) else : x_xz = x[QQ] / (x[QQ] - z[QQ]) alpha = x_xz.min() #end x += alpha * (z - x) iabs = find(abs(x) < tol) ip = find(P[iabs] != -1) ij = iabs[ip] Z[ij] = numpy.array(ij) P[ij] = -1 PP = find(P != -1) ZZ = find(Z != -1) if len(PP) == 1 : XtyPP = Xty[PP] XtXPP = XtX[PP, PP] z[PP] = XtyPP / XtXPP else : XtyPP = numpy.array(Xty[PP]) XtXPP = numpy.array(XtX[PP, numpy.array(PP)[:, numpy.NewAxis]]) z[PP] = numpy.dot(XtyPP, la.generalized_inverse(XtXPP)) #endif z[ZZ] = 0 #end while x = numpy.array(z) w = Xty - numpy.dot(XtX, x) #end while return x, w #end def if __name__ == '__main__' : # # test [x, w] = fnnls(Xt.X, Xt.y, tol) # to solve min ||y - X.x|| s.t. x >= 0 # # matlab:lsqnonneg # X = [1, 10, 4, 10; 4, 5, 1, 12; 5, 1, 9, 20]; # y = [4; 7; 4] # x = lsqnonneg(X, y) => x = [0.9312; 0.3683; 0; 0]; # X = numpy.array([[1, 10, 4, 10], [4, 5, 1, 12], [5, 1, 9, 20]], numpy.Float) y = numpy.array([4, 7, 4], numpy.Float) if False : X = numpy.array([[1, 10, 4, 10], [4, 5, 1, 12], [5, 1, 9, 20], [4, 3, 2, 1]], numpy.Float) y = numpy.array([4, 7, 4, 1], numpy.Float) #end if False : X = numpy.zeros((20, 20), numpy.Float) for n in range(20) : X[n,:] = numpy.arange(0.0, 400.0, step = 20) y = numpy.arange(0.0, 20.0) #end Xt = numpy.transpose(numpy.array(X)) x, w = fnnls(numpy.dot(Xt, X), numpy.dot(Xt, y)) print 'X = ', X print 'y = ', y print 'x = ', x #end if __name__ == '__main__' --Apple-Mail-27-613260171 Content-Transfer-Encoding: quoted-printable Content-Type: text/html; charset=US-ASCII

I've used this to linearise = problems such as:
fitting Neutron TOF data with a = E0 energy spectrum of gaussians (with energy-resolutions a function of = TOF =3D> E): exp(-((E - E0) / deltaE)^2)
fitting = Cerebral Blood Flow data with a k-spectrum set of exponentials: = exp(-kt)
fitting BGO detector pulses (to = detect/correct pile-up) with Tdelay-spectrum 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)-9496-5767
Fax: (613) 9458-5023
=

= --Apple-Mail-27-613260171--