[Scipy-svn] r4754 - in trunk/scipy/optimize: . nnls tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Sep 29 16:15:37 EDT 2008
Author: uwe.schmitt
Date: 2008-09-29 15:13:24 -0500 (Mon, 29 Sep 2008)
New Revision: 4754
Added:
trunk/scipy/optimize/nnls.py
trunk/scipy/optimize/nnls/
trunk/scipy/optimize/nnls/NNLS.F
trunk/scipy/optimize/nnls/nnls.pyf
trunk/scipy/optimize/tests/test_nnls.py
Modified:
trunk/scipy/optimize/SConscript
trunk/scipy/optimize/setup.py
Log:
added nnls, incl. test script
Modified: trunk/scipy/optimize/SConscript
===================================================================
--- trunk/scipy/optimize/SConscript 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/SConscript 2008-09-29 20:13:24 UTC (rev 4754)
@@ -73,6 +73,10 @@
src = [pjoin('minpack2', i) for i in ['dcsrch.f', 'dcstep.f', 'minpack2.pyf']]
env.NumpyPythonExtension('minpack2', source = src)
+# _nnls pyextension
+src = [pjoin('nnls', i) for i in ['NNLS.f', 'nnls.pyf']]
+env.NumpyPythonExtension('_nnls', source = src)
+
# moduleTNC pyextension
env.NumpyPythonExtension('moduleTNC',
source = [pjoin('tnc', i) for i in \
Added: trunk/scipy/optimize/nnls/NNLS.F
===================================================================
--- trunk/scipy/optimize/nnls/NNLS.F 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/nnls/NNLS.F 2008-09-29 20:13:24 UTC (rev 4754)
@@ -0,0 +1,477 @@
+C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
+C
+C Algorithm NNLS: NONNEGATIVE LEAST SQUARES
+C
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 15, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+c
+C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN
+C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
+C
+C A * X = B SUBJECT TO X .GE. 0
+C ------------------------------------------------------------------
+c Subroutine Arguments
+c
+C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE
+C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N
+C MATRIX, A. ON EXIT A() CONTAINS
+C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN
+C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY
+C THIS SUBROUTINE.
+C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON-
+C TAINS Q*B.
+C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL
+C CONTAIN THE SOLUTION VECTOR.
+C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE
+C RESIDUAL VECTOR.
+C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN
+C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0.
+C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z
+C ZZ() AN M-ARRAY OF WORKING SPACE.
+C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N.
+C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS
+C P AND Z AS FOLLOWS..
+C
+C INDEX(1) THRU INDEX(NSETP) = SET P.
+C INDEX(IZ1) THRU INDEX(IZ2) = SET Z.
+C IZ1 = NSETP + 1 = NPP1
+C IZ2 = N
+C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING
+C MEANINGS.
+C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY.
+C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD.
+C EITHER M .LE. 0 OR N .LE. 0.
+C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS.
+C
+C ------------------------------------------------------------------
+ SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE)
+C ------------------------------------------------------------------
+ integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L
+ integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY
+c integer INDEX(N)
+c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M)
+ integer INDEX(*)
+ double precision A(MDA,*), B(*), W(*), X(*), ZZ(*)
+ double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM
+ double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX
+ double precision ZERO, ZTEST
+ parameter(FACTOR = 0.01d0)
+ parameter(TWO = 2.0d0, ZERO = 0.0d0)
+C ------------------------------------------------------------------
+ MODE=1
+ IF (M .le. 0 .or. N .le. 0) then
+ MODE=2
+ RETURN
+ endif
+ ITER=0
+ ITMAX=3*N
+C
+C INITIALIZE THE ARRAYS INDEX() AND X().
+C
+ DO 20 I=1,N
+ X(I)=ZERO
+ 20 INDEX(I)=I
+C
+ IZ2=N
+ IZ1=1
+ NSETP=0
+ NPP1=1
+C ****** MAIN LOOP BEGINS HERE ******
+ 30 CONTINUE
+C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
+C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
+C
+ IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350
+C
+C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
+C
+ DO 50 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ SM=ZERO
+ DO 40 L=NPP1,M
+ 40 SM=SM+A(L,J)*B(L)
+ W(J)=SM
+ 50 continue
+C FIND LARGEST POSITIVE W(J).
+ 60 continue
+ WMAX=ZERO
+ DO 70 IZ=IZ1,IZ2
+ J=INDEX(IZ)
+ IF (W(J) .gt. WMAX) then
+ WMAX=W(J)
+ IZMAX=IZ
+ endif
+ 70 CONTINUE
+C
+C IF WMAX .LE. 0. GO TO TERMINATION.
+C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
+C
+ IF (WMAX .le. ZERO) go to 350
+ IZ=IZMAX
+ J=INDEX(IZ)
+C
+C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
+C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
+C NEAR LINEAR DEPENDENCE.
+C
+ ASAVE=A(NPP1,J)
+ CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0)
+ UNORM=ZERO
+ IF (NSETP .ne. 0) then
+ DO 90 L=1,NSETP
+ 90 UNORM=UNORM+A(L,J)**2
+ endif
+ UNORM=sqrt(UNORM)
+ IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then
+C
+C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ
+C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
+C
+ DO 120 L=1,M
+ 120 ZZ(L)=B(L)
+ CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1)
+ ZTEST=ZZ(NPP1)/A(NPP1,J)
+C
+C SEE IF ZTEST IS POSITIVE
+C
+ IF (ZTEST .gt. ZERO) go to 140
+ endif
+C
+C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
+C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
+C COEFFS AGAIN.
+C
+ A(NPP1,J)=ASAVE
+ W(J)=ZERO
+ GO TO 60
+C
+C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM
+C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER
+C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN
+C COL J, SET W(J)=0.
+C
+ 140 continue
+ DO 150 L=1,M
+ 150 B(L)=ZZ(L)
+C
+ INDEX(IZ)=INDEX(IZ1)
+ INDEX(IZ1)=J
+ IZ1=IZ1+1
+ NSETP=NPP1
+ NPP1=NPP1+1
+C
+ IF (IZ1 .le. IZ2) then
+ DO 160 JZ=IZ1,IZ2
+ JJ=INDEX(JZ)
+ CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1)
+ 160 continue
+ endif
+C
+ IF (NSETP .ne. M) then
+ DO 180 L=NPP1,M
+ 180 A(L,J)=ZERO
+ endif
+C
+ W(J)=ZERO
+C SOLVE THE TRIANGULAR SYSTEM.
+C STORE THE SOLUTION TEMPORARILY IN ZZ().
+ RTNKEY = 1
+ GO TO 400
+ 200 CONTINUE
+C
+C ****** SECONDARY LOOP BEGINS HERE ******
+C
+C ITERATION COUNTER.
+C
+ 210 continue
+ ITER=ITER+1
+ IF (ITER .gt. ITMAX) then
+ MODE=3
+ write (*,'(/a)') ' NNLS quitting on iteration count.'
+ GO TO 350
+ endif
+C
+C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
+C IF NOT COMPUTE ALPHA.
+C
+ ALPHA=TWO
+ DO 240 IP=1,NSETP
+ L=INDEX(IP)
+ IF (ZZ(IP) .le. ZERO) then
+ T=-X(L)/(ZZ(IP)-X(L))
+ IF (ALPHA .gt. T) then
+ ALPHA=T
+ JJ=IP
+ endif
+ endif
+ 240 CONTINUE
+C
+C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
+C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
+C
+ IF (ALPHA.EQ.TWO) GO TO 330
+C
+C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO
+C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
+C
+ DO 250 IP=1,NSETP
+ L=INDEX(IP)
+ X(L)=X(L)+ALPHA*(ZZ(IP)-X(L))
+ 250 continue
+C
+C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
+C FROM SET P TO SET Z.
+C
+ I=INDEX(JJ)
+ 260 continue
+ X(I)=ZERO
+C
+ IF (JJ .ne. NSETP) then
+ JJ=JJ+1
+ DO 280 J=JJ,NSETP
+ II=INDEX(J)
+ INDEX(J-1)=II
+ CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II))
+ A(J,II)=ZERO
+ DO 270 L=1,N
+ IF (L.NE.II) then
+c
+c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
+c
+ TEMP = A(J-1,L)
+ A(J-1,L) = CC*TEMP + SS*A(J,L)
+ A(J,L) =-SS*TEMP + CC*A(J,L)
+ endif
+ 270 CONTINUE
+c
+c Apply procedure G2 (CC,SS,B(J-1),B(J))
+c
+ TEMP = B(J-1)
+ B(J-1) = CC*TEMP + SS*B(J)
+ B(J) =-SS*TEMP + CC*B(J)
+ 280 continue
+ endif
+c
+ NPP1=NSETP
+ NSETP=NSETP-1
+ IZ1=IZ1-1
+ INDEX(IZ1)=I
+C
+C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD
+C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
+C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY
+C THAT ARE NONPOSITIVE WILL BE SET TO ZERO
+C AND MOVED FROM SET P TO SET Z.
+C
+ DO 300 JJ=1,NSETP
+ I=INDEX(JJ)
+ IF (X(I) .le. ZERO) go to 260
+ 300 CONTINUE
+C
+C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK.
+C
+ DO 310 I=1,M
+ 310 ZZ(I)=B(I)
+ RTNKEY = 2
+ GO TO 400
+ 320 CONTINUE
+ GO TO 210
+C ****** END OF SECONDARY LOOP ******
+C
+ 330 continue
+ DO 340 IP=1,NSETP
+ I=INDEX(IP)
+ 340 X(I)=ZZ(IP)
+C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING.
+ GO TO 30
+C
+C ****** END OF MAIN LOOP ******
+C
+C COME TO HERE FOR TERMINATION.
+C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.
+C
+ 350 continue
+ SM=ZERO
+ IF (NPP1 .le. M) then
+ DO 360 I=NPP1,M
+ 360 SM=SM+B(I)**2
+ else
+ DO 380 J=1,N
+ 380 W(J)=ZERO
+ endif
+ RNORM=sqrt(SM)
+ RETURN
+C
+C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE
+C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ().
+C
+ 400 continue
+ DO 430 L=1,NSETP
+ IP=NSETP+1-L
+ IF (L .ne. 1) then
+ DO 410 II=1,IP
+ ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1)
+ 410 continue
+ endif
+ JJ=INDEX(IP)
+ ZZ(IP)=ZZ(IP)/A(IP,JJ)
+ 430 continue
+ go to (200, 320), RTNKEY
+ END
+
+
+ double precision FUNCTION DIFF(X,Y)
+c
+c Function used in tests that depend on machine precision.
+c
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 7, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C
+ double precision X, Y
+ DIFF=X-Y
+ RETURN
+ END
+
+
+C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
+C
+C CONSTRUCTION AND/OR APPLICATION OF A SINGLE
+C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B
+C
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 12, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C ------------------------------------------------------------------
+c Subroutine Arguments
+c
+C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a
+c Householder transformation, or Algorithm H2 to apply a
+c previously constructed transformation.
+C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
+C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO
+C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M
+C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
+C U(),IUE,UP On entry with MODE = 1, U() contains the pivot
+c vector. IUE is the storage increment between elements.
+c On exit when MODE = 1, U() and UP contain quantities
+c defining the vector U of the Householder transformation.
+c on entry with MODE = 2, U() and UP should contain
+c quantities previously computed with MODE = 1. These will
+c not be modified during the entry with MODE = 2.
+C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH
+c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE
+c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED.
+c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS.
+C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
+C ICV STORAGE INCREMENT BETWEEN VECTORS IN C().
+C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
+C NO OPERATIONS WILL BE DONE ON C().
+C ------------------------------------------------------------------
+ SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
+C ------------------------------------------------------------------
+ integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J
+ integer L1, LPIVOT, M, MODE, NCV
+ double precision B, C(*), CL, CLINV, ONE, SM
+c double precision U(IUE,M)
+ double precision U(IUE,*)
+ double precision UP
+ parameter(ONE = 1.0d0)
+C ------------------------------------------------------------------
+ IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
+ CL=abs(U(1,LPIVOT))
+ IF (MODE.EQ.2) GO TO 60
+C ****** CONSTRUCT THE TRANSFORMATION. ******
+ DO 10 J=L1,M
+ 10 CL=MAX(abs(U(1,J)),CL)
+ IF (CL) 130,130,20
+ 20 CLINV=ONE/CL
+ SM=(U(1,LPIVOT)*CLINV)**2
+ DO 30 J=L1,M
+ 30 SM=SM+(U(1,J)*CLINV)**2
+ CL=CL*SQRT(SM)
+ IF (U(1,LPIVOT)) 50,50,40
+ 40 CL=-CL
+ 50 UP=U(1,LPIVOT)-CL
+ U(1,LPIVOT)=CL
+ GO TO 70
+C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ******
+C
+ 60 IF (CL) 130,130,70
+ 70 IF (NCV.LE.0) RETURN
+ B= UP*U(1,LPIVOT)
+C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN.
+C
+ IF (B) 80,130,130
+ 80 B=ONE/B
+ I2=1-ICV+ICE*(LPIVOT-1)
+ INCR=ICE*(L1-LPIVOT)
+ DO 120 J=1,NCV
+ I2=I2+ICV
+ I3=I2+INCR
+ I4=I3
+ SM=C(I2)*UP
+ DO 90 I=L1,M
+ SM=SM+C(I3)*U(1,I)
+ 90 I3=I3+ICE
+ IF (SM) 100,120,100
+ 100 SM=SM*B
+ C(I2)=C(I2)+SM*UP
+ DO 110 I=L1,M
+ C(I4)=C(I4)+SM*U(1,I)
+ 110 I4=I4+ICE
+ 120 CONTINUE
+ 130 RETURN
+ END
+
+
+
+ SUBROUTINE G1 (A,B,CTERM,STERM,SIG)
+c
+C COMPUTE ORTHOGONAL ROTATION MATRIX..
+c
+c The original version of this code was developed by
+c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
+c 1973 JUN 12, and published in the book
+c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
+c Revised FEB 1995 to accompany reprinting of the book by SIAM.
+C
+C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
+C (-S,C) (-S,C)(B) ( 0 )
+C COMPUTE SIG = SQRT(A**2+B**2)
+C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
+C SIG MAY BE IN THE SAME LOCATION AS A OR B .
+C ------------------------------------------------------------------
+ double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO
+ parameter(ONE = 1.0d0, ZERO = 0.0d0)
+C ------------------------------------------------------------------
+ if (abs(A) .gt. abs(B)) then
+ XR=B/A
+ YR=sqrt(ONE+XR**2)
+ CTERM=sign(ONE/YR,A)
+ STERM=CTERM*XR
+ SIG=abs(A)*YR
+ RETURN
+ endif
+
+ if (B .ne. ZERO) then
+ XR=A/B
+ YR=sqrt(ONE+XR**2)
+ STERM=sign(ONE/YR,B)
+ CTERM=STERM*XR
+ SIG=abs(B)*YR
+ RETURN
+ endif
+
+ SIG=ZERO
+ CTERM=ZERO
+ STERM=ONE
+ RETURN
+ END
Added: trunk/scipy/optimize/nnls/nnls.pyf
===================================================================
--- trunk/scipy/optimize/nnls/nnls.pyf 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/nnls/nnls.pyf 2008-09-29 20:13:24 UTC (rev 4754)
@@ -0,0 +1,22 @@
+! -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+python module _nnls ! in
+ interface ! in :_nnls
+ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index_bn,mode) ! in :nnls:NNLS.F
+ double precision dimension(mda,*), intent(copy) :: a
+ integer optional,check(shape(a,0)==mda),depend(a) :: mda=shape(a,0)
+ integer :: m
+ integer :: n
+ double precision dimension(*), intent(copy) :: b
+ double precision dimension(n), intent(out) :: x
+ double precision, intent(out) :: rnorm
+ double precision dimension(*) :: w
+ double precision dimension(*) :: zz
+ integer dimension(*) :: index_bn
+ integer , intent(out):: mode
+ end subroutine nnls
+end python module _nnls
+
+! This file was auto-generated with f2py (version:2_5878).
+! See http://cens.ioc.ee/projects/f2py2e/
Added: trunk/scipy/optimize/nnls.py
===================================================================
--- trunk/scipy/optimize/nnls.py 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/nnls.py 2008-09-29 20:13:24 UTC (rev 4754)
@@ -0,0 +1,41 @@
+import _nnls
+from numpy import asarray_chkfinite, zeros, double
+
+def nnls(A,b, overwrite_a=0, overwrite_b=0):
+ """
+ Solve || Ax - b ||_2 -> min with x>=0
+
+ Inputs:
+ A -- matrix as above
+ b -- vector as above
+
+ Outputs:
+ x -- solution vector
+ rnorm -- residual || Ax-b ||_2
+
+
+ wrapper around NNLS.F code below nnls/ directory
+
+ """
+
+ A,b = map(asarray_chkfinite, (A,b))
+
+ if len(A.shape)!=2:
+ raise ValueError, "expected matrix"
+ if len(b.shape)!=1:
+ raise ValueError, "expected vector"
+
+ m,n = A.shape
+
+ if m != b.shape[0]:
+ raise ValueError, "incompatible dimensions"
+
+ w = zeros((n,), dtype=double)
+ zz = zeros((m,), dtype=double)
+ index=zeros((n,), dtype=int)
+
+ x,rnorm,mode = _nnls.nnls(A,m,n,b,w,zz,index)
+ if mode != 1: raise RuntimeError, "too many iterations"
+
+ return x, rnorm
+
Modified: trunk/scipy/optimize/setup.py
===================================================================
--- trunk/scipy/optimize/setup.py 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/setup.py 2008-09-29 20:13:24 UTC (rev 4754)
@@ -41,6 +41,9 @@
sources = ['slsqp.pyf', 'slsqp_optmz.f']
config.add_extension('_slsqp', sources=[join('slsqp', x) for x in sources])
+ config.add_extension('_nnls', sources=[join('nnls', x) \
+ for x in ["NNLS.f","nnls.pyf"]])
+
config.add_data_dir('tests')
config.add_data_dir('benchmarks')
return config
Added: trunk/scipy/optimize/tests/test_nnls.py
===================================================================
--- trunk/scipy/optimize/tests/test_nnls.py 2008-09-27 21:46:03 UTC (rev 4753)
+++ trunk/scipy/optimize/tests/test_nnls.py 2008-09-29 20:13:24 UTC (rev 4754)
@@ -0,0 +1,31 @@
+""" Unit tests for nonlinear solvers
+Author: Ondrej Certik
+May 2007
+"""
+
+from numpy.testing import *
+
+from scipy.optimize import nnls
+from numpy import arange, dot
+from numpy.linalg import norm
+
+
+class TestNNLS(TestCase):
+ """ Test case for a simple constrained entropy maximization problem
+ (the machine translation example of Berger et al in
+ Computational Linguistics, vol 22, num 1, pp 39--72, 1996.)
+ """
+
+ def test_nnls(self):
+ a=arange(25.0).reshape(-1,5)
+ x=arange(5.0)
+ y=dot(a,x)
+ x, res= nnls.nnls(a,y)
+ assert res<1e-7
+ assert norm(dot(a,x)-y)<1e-7
+
+if __name__ == "__main__":
+ run_module_suite()
+
+
+
More information about the Scipy-svn
mailing list