[Scipy-svn] r4754 - in trunk/scipy/optimize: . nnls tests

scipy-svn@scip... scipy-svn@scip...
Mon Sep 29 15:15:37 CDT 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