Skip to content

Loading…

ENH: Add bounded linear least squares routine. Closes #3130. #3137

Closed
wants to merge 17 commits into from
@jseabold

This works, but I'd like to get some input from f2py experts (cc // @pv) on a few points. You can run the code with the below. The Fortran code is hard-coded for float32 precision. Presumably, that shouldn't be a big deal to change. However, I'm running into some other issues that are I think related to using Fortran 90. I've never wrapped Fortran 90 code before.

I can't change the intent of, say, a to be (in, out, copy). If I do, I receive an error like this from f2py. That's the main sticking point for me, I think. However, the wrapping code makes a copy, so maybe we can live with it?

gfortran:f77: /tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f
/tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f:27.31:

          real(kind=kind(one)), dimension(:,:),intent(in,out,copy)      
                               1
Error: Invalid character in name at (1)
/tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f:39.16:

      call bvls(a, b, bnd, x, rnorm, nsetp, w, index_bn, ierr)          
                1
Error: Rank mismatch in argument 'a' at (1) (scalar and rank-2)
/tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f:27.31:

          real(kind=kind(one)), dimension(:,:),intent(in,out,copy)      
                               1
Error: Invalid character in name at (1)
/tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f:39.16:

      call bvls(a, b, bnd, x, rnorm, nsetp, w, index_bn, ierr)          
                1
Error: Rank mismatch in argument 'a' at (1) (scalar and rank-2)
error: Command "/usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops -I/tmp/tmp_TCHXo/src.linux-x86_64-2.7 -I/usr/local/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c -c /tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.f -o /tmp/tmp_TCHXo/tmp/tmp_TCHXo/src.linux-x86_64-2.7/_bvls-f2pywrappers.o" failed with exit status 1
import numpy as np
from scipy import optimize

x = np.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
                  -2.8394297E-01,1.3416969E+00,1.3757038E+00,
                  -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
                  8.2065094E-01])[:,None]
x = np.column_stack((np.ones_like(x), x))
np.random.seed(12)
y = np.dot(x, [1.8, 5.4]) + np.random.random(len(x))
bounds = None

print optimize.bounded_lstsq(x, y, bounds=bounds)

bounds = [(None, None), (None, 4)]
print optimize.bounded_lstsq(x, y, bounds=bounds)
@josef-pkt
SciPy member

I thought the issue says linear least squares, lstsq not leastsq

@jseabold

It does. I didn't know this was a convention. Can change it before its merged.

@jseabold

Typo in the commit and PR...

@coveralls

Coverage Status

Coverage remained the same when pulling fa9258c7d9316471bd1960f473ea635f4ab4674b on jseabold:constrained-lstsq into d84e489 on scipy:master.

@josef-pkt
SciPy member

linalg.lstsq and optimize.leastsq are just the current names, not much of a convention
but using leastsq in the name for a linear least squares would be confusing to me.

maybe bvls in analogy to nnls would be a better name, or something more verbose

(I've no idea about the Fortran, obviously.)

@coveralls

Coverage Status

Coverage remained the same when pulling fa9258c7d9316471bd1960f473ea635f4ab4674b on jseabold:constrained-lstsq into d84e489 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling d5e59db on jseabold:constrained-lstsq into d84e489 on scipy:master.

@pv pv commented on an outdated diff
scipy/optimize/bvls.py
((45 lines not shown))
+
+ if A.ndim != 2:
+ raise ValueError("A must be 2-dimensional")
+ m, n = A.shape
+
+ if b.ndim != 1:
+ raise VelueError("Expected 1-dimensional b")
+ if A.shape[0] != b.shape[0]:
+ raise ValueError("A and b are not conformable")
+
+ if bounds is None or len(bounds) == 0:
+ bnds = array([[-1.0E12]*n, [1.0E12]*n])
+ else:
+ bnds = array(bounds, float).T
+ infbnd = ~isfinite(bnds)
+ bnds[0, infbnd[0]] = -inf
@pv SciPy member
pv added a note

Not needed, I believe.

@jseabold
jseabold added a note

It allows using None for bounds. Passing NaNs as Infs doesn't work I don't think.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on an outdated diff
scipy/optimize/bvls.py
((41 lines not shown))
+ A, b = map(asarray_chkfinite, (A, b))
+ #TODO: remove this when sorted
+ A = A.astype(dtype=float32, order='F')
+ b = b.astype(dtype=float32, order='F')
+
+ if A.ndim != 2:
+ raise ValueError("A must be 2-dimensional")
+ m, n = A.shape
+
+ if b.ndim != 1:
+ raise VelueError("Expected 1-dimensional b")
+ if A.shape[0] != b.shape[0]:
+ raise ValueError("A and b are not conformable")
+
+ if bounds is None or len(bounds) == 0:
+ bnds = array([[-1.0E12]*n, [1.0E12]*n])
@pv SciPy member
pv added a note

Why not -inf, inf?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on an outdated diff
scipy/optimize/bvls.py
((6 lines not shown))
+def bounded_lstsq(A, b, bounds=()):
+ """
+ Bounded variable (linear) least-squares.
+
+ Solve ``argmin_x || Ax - b ||_2`` for possibly bounded ``x``.
+
+ Parameters
+ ----------
+ A : ndarray
+ Matrix ``A`` as shown above.
+ b : ndarray
+ Right-hand side vector.
+ bounds : list
+ A list of tuples specifying the lower and upper bound for each
+ independent variables [(xl0, xu0), (xl1, xu1), ...] Infinite values
+ or None will be interpreted as large floating values, negated as
@pv SciPy member
pv added a note

This is not what the code actually does?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on an outdated diff
scipy/optimize/bvls/bvls.f90
@@ -0,0 +1,703 @@
+SUBROUTINE BVLS ( A, B, BND, X, RNORM, NSETP, W, INDEX, IERR )
@pv SciPy member
pv added a note

The source file should state that we have obtained a license from the authors to use it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on an outdated diff
scipy/optimize/bvls.py
@@ -0,0 +1,79 @@
+from . import _bvls
+from numpy import asarray_chkfinite, zeros, float32, array, isfinite, inf
+
+__all__ = ['bounded_lstsq']
+
+def bounded_lstsq(A, b, bounds=()):
@pv SciPy member
pv added a note

Or, lstsq_bounded?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv
SciPy member
pv commented

The f2py part looks OK to me.

What on the other hand doesn't look so nice is that it uses float32 rather than float64. This could probably be changed in the source.

@keflavich

Any chance this will be merged soon? I've wanted this feature in scipy for years!

@jseabold

It needs a bit more work. Will get to it at some point when I have an afternoon, feel free to submit a new PR if you beat me to it.

@keflavich keflavich referenced this pull request in jseabold/scipy
Closed

Attempt to fix up issues in bvls #2

@jseabold

Converted float32 -> float64 in the fortran. Tested unbounded vs. numpy.linalg.lstsq. Addressed @pv's comments except I didn't change the name, but I will change it as desired. Rebased on master.

@coveralls

Coverage Status

Coverage remained the same when pulling c447034 on jseabold:constrained-lstsq into 9ccc684 on scipy:master.

@coveralls

Coverage Status

Coverage remained the same when pulling c447034 on jseabold:constrained-lstsq into 9ccc684 on scipy:master.

@jseabold

The remaining issue is for numpy 1.5.1. What's the recommended way to do

A.astype(float, order='F')

in 1.5.1? This?

A.astype(float).copy('F')

Is that two copies now...?

@jseabold jseabold closed this
@jseabold jseabold reopened this
@coveralls

Coverage Status

Coverage remained the same when pulling e8b0631 on jseabold:constrained-lstsq into cb946a5 on scipy:master.

@jseabold

Ok, I think I finally wrote code that should work on numpy 1.5.1 as well. It's returning an error code that indicates that the shape of A is not 2-d. I have no idea how this can be. Have there been any changes to f2py since numpy 1.5.1 (last 3 years) that might lead to this behavior? I have no idea.

If I can pin down the last numpy release this worked on, can I 1) make the test conditional or 2) make the function conditionally defined?

@rgommers
SciPy member

Not many changes in f2py, a handful of bug fixes and a few new features. For 1.5.1, does it not work at all or does it fail with some inputs?

@rgommers
SciPy member

Why would we be able to merge Fortran90 code here? If it doesn't work in scipy.integrate, it shouldn't work here.....

@jseabold

Like I said, I don't know anything about Fortran 90 and f2py. I'd be interested in some background here. The code is simple enough that I don't think removing those language features that aren't ok would be a big deal. It's just mysterious to me that it fails for numpy 1.5.1. I tried to do a git bisect last night but I failed. Virtualenvs / I don't know how to nosetest within source of scipy, apparently. Long story.

Somehow one or both of size(A, 1) of size(A, 2) return <= 0 on numpy 1.5.1, for any input as far as I can tell, though I don't have time to build another old numpy right now.

@WarrenWeckesser

Fortran 90/95 issue: #2829

@WarrenWeckesser WarrenWeckesser commented on the diff
scipy/optimize/bento.info
@@ -41,3 +41,7 @@ Library:
Sources:
slsqp/slsqp_optmz.f,
slsqp/slsqp.pyf
+ Extension: _bvls
+ Sources:
+ bvls/bvls.f,

bvls.f90

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@WarrenWeckesser

Unfortunately, REAL(KIND=8) is not guaranteed to mean DOUBLE PRECISION for all compilers. The interpretation of the <int> in KIND=<int> is compiler dependent.

A portable method is to use SELECTED_REAL_KIND(...) (e.g. http://fortranwiki.org/fortran/show/selected_real_kind), or do something like the suggestion here: http://www.fortran90.org/src/best-practices.html#floating-point-numbers

See also http://stackoverflow.com/questions/12523524/f2py-specifying-real-precision-in-fortran-when-interfacing-with-python

@jseabold

Ok, thanks @WarrenWeckesser for the pointers. Fun stuff. It sounds like selected_real_kind is a no-go anyway because we need to stick with F77 because of 64-bit windows.

F2py should work fine with some fairly broad subset of f90. I went through and removed most of the f90 issues during the first round of cleanup I think. Though obviously not completely. So the fortran 90 issues on windows and compilers issues shouldn't explain the current failure for only 1.5.1 on travis.

I'll need to clean up the code anyway for F77, so I'll see if it falls out doing this.

@jseabold

I guess an implicit question in the above is is it worth translating f90 -> f77 or waiting for cross-platform f90 compilation? Translating seems kind of ridiculous when I step back for a second.

@jseabold

That seems to have solved the 1.5.1 issues. The question now is how do I check if the code is fully f77 compliant? Current gfortran only allows automatic checking back to f95 it looks like. Any ideas?

Then I'll have to work on the cross-platorm/compiler type declarations

@coveralls

Coverage Status

Coverage remained the same when pulling 56dbaa9 on jseabold:constrained-lstsq into cb946a5 on scipy:master.

@rgommers
SciPy member

@jseabold we can just check if it works with the scipy release build setup, which uses g77.

@jseabold

The latest commit should make it so that double precision is used consistently across platforms and compilers.

@jseabold

My last thought was whether this makes sense in optimize or if it should go into scipy.linalg to mirror numpy.linalg.lstsq?

@rgommers
SciPy member

+1 for optimize

@jseabold

FYI, the builds are fine the full test suite just timed out. Let me know how to go about testing with g77 or what the outcomes of this are. I suspect some more changes may need to be made.

@Tillsten

Because nnls is already in optimize, it should also go into it.

@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((28 lines not shown))
+ rnorm : float
+ The residual, ``|| Ax-b ||_2``.
+
+ See also
+ --------
+ scipy.optimize.nnls
+
+ Notes
+ -----
+ The Fortran code for this was originally written by Charles Lawson and
+ Richard Hanson, who agreed to release the code under the BSD license
+ for inclusion in scipy.
+
+ Examples
+ --------
+ >>> import numpy as np
@rgommers SciPy member

this import can be left out, standard convention

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((43 lines not shown))
+ >>> import numpy as np
+ >>> from scipy import optimize
+
+ >>> x = np.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
+ ... -2.8394297E-01,1.3416969E+00,1.3757038E+00,
+ ... -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
+ ... 8.2065094E-01])[:,None]
+ >>> x = np.column_stack((np.ones_like(x), x))
+ >>> np.random.seed(12)
+ >>> y = np.dot(x, [1.8, 5.4]) + np.random.random(len(x))
+ >>> bounds = None
+
+ >>> print optimize.bounded_lstsq(x, y, bounds=bounds)
+
+ >>> bounds = [(None, None), (None, 4)]
+ >>> print optimize.bounded_lstsq(x, y, bounds=bounds)
@rgommers SciPy member

Please no print in docstrings (also above). This way is invalid syntax in Python 3.x, and prints aren't needed in an interpreter.

Also missing the actual output.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
@@ -0,0 +1,98 @@
+from . import _bvls
+from numpy import zeros, array, isfinite, inf
@rgommers SciPy member

This is new code, would be nicer to not do this and use the np prefix.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((6 lines not shown))
+
+def bounded_lstsq(A, b, bounds=()):
+ """
+ Bounded variable (linear) least-squares.
+
+ Solve ``argmin_x || Ax - b ||_2`` for possibly bounded ``x``.
+
+ Parameters
+ ----------
+ A : ndarray
+ Matrix ``A`` as shown above.
+ b : ndarray
+ Right-hand side vector.
+ bounds : list
+ A list of tuples specifying the lower and upper bound for each
+ independent variables [(xl0, xu0), (xl1, xu1), ...]. None or
@rgommers SciPy member

Need double backticks for this expression.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/__init__.py
@@ -172,6 +172,7 @@
from .nonlin import *
from .slsqp import fmin_slsqp
from .nnls import nnls
+from .bvls import bounded_lstsq
@rgommers SciPy member

bvls.py should be _bvls.py

@rgommers SciPy member

and then the compiled extension called something else. Or rename the Python file to something else, but needs to be private.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((4 lines not shown))
+__all__ = ['bounded_lstsq']
+
+
+def bounded_lstsq(A, b, bounds=()):
+ """
+ Bounded variable (linear) least-squares.
+
+ Solve ``argmin_x || Ax - b ||_2`` for possibly bounded ``x``.
+
+ Parameters
+ ----------
+ A : ndarray
+ Matrix ``A`` as shown above.
+ b : ndarray
+ Right-hand side vector.
+ bounds : list
@rgommers SciPy member

You can also allow an array of shape (N, 2). Already works with this code.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((62 lines not shown))
+ if not isfinite(A).all() or not isfinite(b).all():
+ raise ValueError("A and b may not contain NaNs or infs")
+
+ if A.ndim != 2:
+ raise ValueError("A must be 2-dimensional")
+ m, n = A.shape
+
+ if b.ndim != 1:
+ raise ValueError("Expected 1-dimensional b")
+ if A.shape[0] != b.shape[0]:
+ raise ValueError("A and b are not conformable")
+
+ if bounds is None or len(bounds) == 0:
+ bnds = array([[-inf]*n, [inf]*n])
+ else:
+ bnds = array(bounds, float).T
@rgommers SciPy member

If array_like bounds is allowed, this should be np.asarray. And style nitpick, I prefer dtype=float.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((70 lines not shown))
+ raise ValueError("Expected 1-dimensional b")
+ if A.shape[0] != b.shape[0]:
+ raise ValueError("A and b are not conformable")
+
+ if bounds is None or len(bounds) == 0:
+ bnds = array([[-inf]*n, [inf]*n])
+ else:
+ bnds = array(bounds, float).T
+ infbnd = ~isfinite(bnds)
+ bnds[0, infbnd[0]] = -inf
+ bnds[1, infbnd[1]] = inf
+ if bnds.shape[1] != n:
+ raise ValueError("The length of bounds is not compatible with "
+ "Ax=b. Got %d. Expected %d" (len(bnds), n))
+
+ w = zeros((n,), dtype=float, order='F')
@rgommers SciPy member

zeros(n) instead of zeros((n,))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

@jseabold did you email the original authors for permission? Reading the comments I'm not sure the copyright is theirs or is owned by SIAM. And IIRC that Burkardt site isn't very reliable with licensing of redistributed code.

@rkern does this look OK to you?

@rgommers
SciPy member

Did you send a message to the mailing list about this enhancement? If so, I can't find it right now. If not, would be good to do that.

@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((18 lines not shown))
+ Right-hand side vector.
+ bounds : list
+ A list of tuples specifying the lower and upper bound for each
+ independent variables [(xl0, xu0), (xl1, xu1), ...]. None or
+ -inf, inf, can be used to indicate no bounds.
+
+ Returns
+ -------
+ x : ndarray
+ Solution vector
+ rnorm : float
+ The residual, ``|| Ax-b ||_2``.
+
+ See also
+ --------
+ scipy.optimize.nnls
@rgommers SciPy member

remove scipy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((24 lines not shown))
+ Returns
+ -------
+ x : ndarray
+ Solution vector
+ rnorm : float
+ The residual, ``|| Ax-b ||_2``.
+
+ See also
+ --------
+ scipy.optimize.nnls
+
+ Notes
+ -----
+ The Fortran code for this was originally written by Charles Lawson and
+ Richard Hanson, who agreed to release the code under the BSD license
+ for inclusion in scipy.
@rgommers SciPy member

I'd suggest to put this in a comment at the top of the file, and not in the docstring.

@rgommers SciPy member

This section could benefit from some notes on the algorithm or a link to a paper.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/bvls.py
((82 lines not shown))
+ raise ValueError("The length of bounds is not compatible with "
+ "Ax=b. Got %d. Expected %d" (len(bnds), n))
+
+ w = zeros((n,), dtype=float, order='F')
+ index = zeros((n,), dtype=int, order='F')
+ x = zeros((n,), dtype=float, order='F')
+
+ rnorm, nsetp, ierr = _bvls.bvls(A, b, bnds, x, w, index)
+ if ierr == 1:
+ raise ValueError("M <= 0 or N <= 0")
+ elif ierr == 2:
+ raise ValueError("Size or shape violation.")
+ elif ierr == 3:
+ raise ValueError("Input bounds are inconsistent")
+ elif ierr == 4:
+ raise ValueError("Exceeded maximum number of iterations.")
@rgommers SciPy member

A dict of error message strings, and a single if ierr > 0: raise ValueError(_errmsg[ierr]) would be cleaner. That's how it's done in a number of other places in scipy.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on an outdated diff
scipy/optimize/tests/test_bvls.py
@@ -0,0 +1,30 @@
+from scipy.optimize import bounded_lstsq
+from numpy import arange, dot, random
+from numpy.linalg import norm, lstsq
@rgommers SciPy member

Same comment as earlier, I'd prefer import numpy as np.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

You need to rename bvls.f90 for g77 otherwise it just refuses to try. After doing that, the MinGW build with g77 fails: https://gist.github.com/rgommers/9590685

@rkern
SciPy member

@rgommers I believe that permission was obtained in 2008 (see the attachment).

@rgommers
SciPy member

@rkern thanks!

@jseabold

Thanks for the comments, Ralf. I addressed them all. It builds locally for me using .f extension, but I don't know how strict the compatibility checks are. The tests also pass locally. I'll ping the ML. Any specific concerns or just to ask for comments?

@Tillsten

I would really like to see this going into 0.14, any chance?

@rgommers
SciPy member

@Tillsten no chance I'm afraid. It's hard enough getting reliable releases out the door without adding new functionality (and especially C/Fortran extensions) after the beta.

@rgommers
SciPy member

@jseabold thanks, changes all look good. I'll test again with MinGW / g77.

@WarrenWeckesser

There is still some more work needed to get bvls.f to compile with g77.

Declarations using the :: syntax are not part of Fortran 77. A line like this

  real*8, parameter :: ZERO = 0E0, ONE=1E0, TWO = 2E0

should be changed to something like

  real*8 ZERO, ONE, TWO
  parameter (ZERO=0D0, ONE=1D0, TWO=2D0)

(Also note the change from E to D, indicating double precision constants.)

Fortran 77 doesn't have assumed-shape arrays, so a declaration such as

     REAL*8 U(:)

will not compile.

CONTAINS is part of the Fortran 90 module interface specification language, so this line:

  CONTAINS  ! These are internal subroutines.

must be removed.

I'm using a Mac, and I just tried compiling it with the g77 compiler from http://hpc.sourceforge.net/. Those were just a few of the problems reported.

It looks like there are less than 400 lines of non-comment Fortran code. Perhaps, instead of incrementally converting to Fortran 77 until it compiles with g77 (a sketchy strategy to begin with), it would make more sense to invest the time in converting it to C or Cython? (That's assuming we don't expect to have a Fortran 90 compiler for Windows any time soon.)

@rgommers
SciPy member

There's not much movement on that front, so I wouldn't count on being able to drop g77 any time soon:(

@jseabold

I agree that the incremental update approach is probably less than ideal. The bonuses are I don't have to think about it, and it's less time-consuming for me (I can do it in the background while I'm working and fix it until it isn't broken). I'm not going to translate this to Cython anytime soon. Maybe f2c would be appropriate here? I've never used it, but I'm generally more comfortable in C than Fortran (77).

Alternatively, can you post the build command you used with g77, so I don't have to rely on y'all running and reporting back? I'll see if I can find an old g77 for my machine. From what I understand the g77 that ships with f2c on linux is not the same thing as old gcc g77. Though I don't know for sure.

@WarrenWeckesser

Alternatively, can you post the build command you used with g77

I used just the basic command g77 -c bvls.f.

@sturlamolden

Why bring more Fortran into existance? If it's for the sake of LAPACK, use LAPACKE from C or Cython.

@sturlamolden

f2c and g77 was written by the same guy. It is basically the same compiler except for the backend. There are also two scripts called f77 and fort77 that does the same as g77 except with f2c and cc. If it works on f2c you can be pretty sure it works on g77. Today a modern C compiler will probably optimize better than g77, so we should look into using f2c instead on Windows.

@sturlamolden sturlamolden commented on the diff
scipy/optimize/bvls/bvls.f
((514 lines not shown))
+ IBOUND=1
+ EXIT
+ ELSE IF ( X(I).ge.BND(2,I)) then
+ IBOUND=2
+ EXIT
+ END IF
+ END DO
+ IF (IBOUND.le.0) EXIT
+ END DO
+!
+! Copy B( ) into Z( ). Then solve again and loop back.
+!
+ Z(1:M)=B(1:M)
+
+ DO I = NSETP, 1, -1
+ IF (I.ne.NSETP) Z(1:I)=Z(1:I)-A(1:I,II)*Z(I+1)

Vectorized expressions are invalid in Fortran 77.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@sturlamolden sturlamolden commented on the diff
scipy/optimize/bvls/bvls.f
((558 lines not shown))
+ L=INDEX(IP)
+ IF (Z(IP).le.BND(1,L)) then
+!
+! Z(IP) HITS LOWER BOUND
+!
+ LBOUND=1
+ ELSE IF (Z(IP).ge.BND(2,L)) then
+!
+! Z(IP) HITS UPPER BOUND
+!
+ LBOUND=2
+ else
+ LBOUND = 0
+ END IF
+
+ IF ( LBOUND.ne.0 ) then

Using name of intrinsic function as variable name.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@sturlamolden sturlamolden commented on the diff
scipy/optimize/bvls/bvls.f
((325 lines not shown))
+! Set FREE2 = true if X(J) is not at the right end-point of its
+! constraint region.
+! Set FREE = true if X(J) is not at either end-point of its
+! constraint region.
+!
+ FREE1 = X(J).gt.BND(1,J)
+ FREE2 = X(J).lt.BND(2,J)
+ FREE = FREE1 .and. FREE2
+
+ IF ( FREE ) then
+ CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+ else
+!
+! Compute dual coefficient W(J).
+!
+ W(J)=dot_product(A(NPP1:M,J),B(NPP1:M))

Use BLAS *DOT instead of intrinsic fortran 90 dot_product

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@sturlamolden

The Fortran code is littered with Fortran 90 syntax and intrinsic functions. Transforming to Fortran 77 cannot be done easily. Also the Fortran 90 code has errors like hard-coded kind numbers (allowed in FORTRAN IV, deprecated in Fortran 77). I vote for not including this code in SciPy. The Fortran is shitty and needs a major overhaul. Then it could just as well be rewritten in C or Cython.

@jseabold

Ok, this PR is dead until I (someone) rewrites this in C/Cython. I don't have any time that I can earmark for this anytime soon I don't think.

@jseabold

Thanks @sturlamolden for the sanity check.

@newville

Just curious -- was dqed.f from netlib.org/opt (which is Fortran77) considered? I've no experience with this code, but it seems to solve essentially the same problem.

@jseabold

That's news to me. I'll look into it. (And if that's what the Fortran 77 version of this code would look like, then yikes.)

@Tillsten

deqed is not optimized to solve linear least squares.

@jseabold

Meaning that it does it, but only as a special case so it's inefficient?

@Tillsten

Yep, bvls and most other bounded linear least squares solvers - as far as i can tell - are using an active set method, where the free sub problem is solved via linear least squares. This should be faster for most cases.

As already mentioned in the issue, https://bitbucket.org/tillsten/pymls contains a python implementation of a matlab package. But it is neither tested well nor did i fully understand what i was implementing. I think someone with a better understanding of the method could do a good python implementation in a day.

@pv pv removed the PR label
@ev-br
SciPy member

gh-5110 adds a python implementation of BVLS.
cc @nmayorov

@jseabold
@charris
SciPy member

@jseabold Will do. Thanks for the work.

@charris charris closed this
@ev-br
SciPy member

gh-5110 has landed, would be great to kick the tires

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
1 .gitignore
@@ -146,6 +146,7 @@ scipy/optimize/cobyla/_cobylamodule.c
scipy/optimize/lbfgsb/_lbfgsbmodule.c
scipy/optimize/minpack2/minpack2module.c
scipy/optimize/nnls/_nnlsmodule.c
+scipy/optimize/bvls/_bvlsmodule.c
scipy/optimize/slsqp/_slsqpmodule.c
scipy/signal/_spectral.c
scipy/signal/correlate_nd.c
View
1 scipy/linalg/basic.py
@@ -500,6 +500,7 @@ def lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False,
See Also
--------
optimize.nnls : linear least squares with non-negativity constraint
+ optimize.bounded_lstsq : linear least squares with arbitrary bounds
"""
View
1 scipy/optimize/__init__.py
@@ -172,6 +172,7 @@
from .nonlin import *
from .slsqp import fmin_slsqp
from .nnls import nnls
+from ._bvls import bounded_lstsq
from ._basinhopping import basinhopping
__all__ = [s for s in dir() if not s.startswith('_')]
View
115 scipy/optimize/_bvls.py
@@ -0,0 +1,115 @@
+"""
+The Fortran code for this was originally written by Charles Lawson and
+Richard Hanson, who agreed to release the code under the BSD license
+for inclusion in scipy.
+"""
+
+from . import _bvlslib
+import numpy as np
+
+__all__ = ['bounded_lstsq']
+
+_error_msg_dict = {
+ 1: "M <= 0 or N <= 0",
+ 2: "Size or shape violation.",
+ 3: "Input bounds are inconsistent",
+ 4: "Exceeded maximum number of iterations."
+ }
+
+
+def bounded_lstsq(A, b, bounds=()):
+ """
+ Bounded variable (linear) least-squares.
+
+ Solve ``argmin_x || Ax - b ||_2`` for possibly bounded ``x``.
+
+ Parameters
+ ----------
+ A : ndarray
+ Matrix ``A`` as shown above.
+ b : ndarray
+ Right-hand side vector.
+ bounds : array_like
+ An (N,2) shape array or a list of tuples specifying the lower and upper
+ bound for each independent variables ``[(xl0, xu0), (xl1, xu1), ...]``.
+ None or -inf, inf, can be used to indicate no bounds.
+
+ Returns
+ -------
+ x : ndarray
+ Solution vector
+ rnorm : float
+ The residual, ``|| Ax-b ||_2``.
+ nsetp : int
+ The number of items in x that are not at their constraint value
+
+ See also
+ --------
+ scipy.optimize.nnls
+
+ Examples
+ --------
+ >>> from scipy import optimize
+
+ >>> x = np.array([-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
+ ... -2.8394297E-01,1.3416969E+00,1.3757038E+00,
+ ... -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
+ ... 8.2065094E-01])[:,None]
+ >>> x = np.column_stack((np.ones_like(x), x))
+ >>> np.random.seed(12)
+ >>> y = np.dot(x, [1.8, 5.4]) + np.random.random(len(x))
+ >>> bounds = None
+
+ >>> optimize.bounded_lstsq(x, y, bounds=bounds)
+ (array([ 2.26327373, 5.42102807]), 1.160159078139885, 2)
+
+ >>> bounds = [(None, None), (None, 4)]
+ >>> optimize.bounded_lstsq(x, y, bounds=bounds)
+ (array([ 2.39941956, 4. ]), 5.379227084630758, 1)
+
+ References
+ ----------
+
+ The reference for the original Fortran code ::
+
+ Charles Lawson, Richard Hanson,
+ Solving Least Squares Problems,
+ SIAM, 1995,
+ ISBN: 0898713560,
+ LC: QA275.L38.
+ """
+ A = np.array(A, dtype=float, copy=True, order='F')
+ b = np.array(b, dtype=float, copy=True, order='F')
+ if not np.isfinite(A).all() or not np.isfinite(b).all():
+ raise ValueError("A and b may not contain NaNs or infs")
+
+ if A.ndim != 2:
+ raise ValueError("A must be 2-dimensional")
+ m, n = A.shape
+
+ if b.ndim != 1:
+ raise ValueError("Expected 1-dimensional b")
+ if A.shape[0] != b.shape[0]:
+ raise ValueError("A and b are not conformable")
+
+ if bounds is None or len(bounds) == 0:
+ bnds = np.array([[-np.inf]*n, [np.inf]*n])
+ else:
+ bnds = np.asarray(bounds, dtype=float).T
+ infbnd = ~np.isfinite(bnds)
+ bnds[0, infbnd[0]] = -np.inf
+ bnds[1, infbnd[1]] = np.inf
+ if bnds.shape[1] != n:
+ raise ValueError("The length of bounds is not compatible with "
+ "Ax=b. Got %d. Expected %d" (len(bnds), n))
+
+ w = np.zeros(n, dtype=float, order='F')
+ index = np.zeros(n, dtype=int, order='F')
+ x = np.zeros(n, dtype=float, order='F')
+
+ rnorm, nsetp, ierr = _bvlslib.bvls(A, b, bnds, x, w, index)
+
+ if ierr > 0:
+ raise ValueError(_error_msg_dict[ierr])
+
+ return x, rnorm, nsetp
View
4 scipy/optimize/bento.info
@@ -41,3 +41,7 @@ Library:
Sources:
slsqp/slsqp_optmz.f,
slsqp/slsqp.pyf
+ Extension: _bvls
+ Sources:
+ bvls/bvls.f,

bvls.f90

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ bvls/bvls.pyf
View
1 scipy/optimize/bscript
@@ -11,6 +11,7 @@ def pre_build(context):
context.tweak_extension("minpack2", features=f2py_features, use=use)
context.tweak_extension("_nnls", features=f2py_features, use=use)
context.tweak_extension("_slsqp", features=f2py_features, use=use)
+ context.tweak_extension("_bvls", features=f2py_features, use=use)
context.tweak_extension("_zeros", features=features, use=use)
context.tweak_extension("_minpack", features=features, use=use)
View
1 scipy/optimize/bvls/README.md
@@ -0,0 +1 @@
+bvls.f90 was obtained on December 6, 2013 from http://people.sc.fsu.edu/~%20jburkardt/f_src/bvls/bvls.html. It was renamed to bvls.f and modified to build with F77 compilers. The timestamp subroutine, which was under the LGPL license was removed.
View
716 scipy/optimize/bvls/bvls.f
@@ -0,0 +1,716 @@
+ SUBROUTINE BVLS ( A, B, M, N, BND, X, RNORM, NSETP, W, INDEX,
+ + IERR )
+
+!*****************************************************************************80
+!
+!! BVLS solves a least squares problem with bounds on the variables.
+!
+! Given an M by N matrix, A, and an M-vector, B, compute an
+! N-vector, X, that solves the least-squares problem A *X = B
+! subject to X(J) satisfying BND(1,J) <= X(J) <= BND(2,J)
+!
+! The values BND(1,J) = -huge(ONE) and BND(2,J) = huge(ONE) are
+! suggested choices to designate that there is no constraint in that
+! direction. The parameter ONE is 1.0 in the working precision.
+!
+! This algorithm is a generalization of NNLS, that solves
+! the least-squares problem, A * X = B, subject to all X(J) >= 0.
+! The subroutine NNLS appeared in 'SOLVING LEAST SQUARES PROBLEMS,'
+! by Lawson and Hanson, Prentice-Hall, 1974. Work on BVLS was started
+! by C. L. Lawson and R. J. Hanson at Jet Propulsion Laboratory,
+! 1973 June 12. Many modifications were subsequently made.
+! This Fortran 90 code was completed in April, 1995 by R. J. Hanson.
+! Assumed shape arrays, automatic arrays, array operations, internal
+! subroutines and F90 elemental functions are used. Comments are
+! prefixed with ! More than 72 characters appear on some lines.
+! The BVLS package is an additional item for the reprinting of the book
+! by SIAM Publications and the distribution of the code package
+! using netlib and Internet or network facilities.
+!
+! Modified:
+!
+! 20 October 2008
+!
+! Author:
+!
+! Charles Lawson, Richard Hanson
+!
+! Reference:
+!
+! Charles Lawson, Richard Hanson,
+! Solving Least Squares Problems,
+! SIAM, 1995,
+! ISBN: 0898713560,
+! LC: QA275.L38.
+!
+! License:
+!
+! The code was released under the BSD license by the original authors
+! for inclusion into SciPy.
+!
+!
+!INTERFACE
+! SUBROUTINE BVLS (A, B, BND, X, RNORM, NSETP, W, INDEX, IERR)
+! REAL*8 A(:,:), B(:), BND(:,:), X(:), RNORM, W(:)
+! INTEGER NSETP, INDEX(:), IERR
+! END SUBROUTINE
+!END INTERFACE
+!
+! Parameters:
+!
+! A(M,N) [INTENT(InOut)]
+! On entry A() contains the M by N matrix, A.
+! On return A() contains the product matrix, Q*A, where
+! Q is an M by M orthogonal matrix generated by this
+! subroutine. The dimensions are M=size(A,1) and N=size(A,2).
+!
+! B(M) [INTENT(InOut)]
+! On entry B() contains the M-vector, B.
+! On return, B() contains Q*B. The same Q multiplies A.
+!
+! M [INTENT(In)]
+! The leading dimension of A
+!
+! N [INTENT(In)]
+! The second dimension of A
+!
+! BND(2,N) [INTENT(In)]
+! BND(1,J) is the lower bound for X(J).
+! BND(2,J) is the upper bound for X(J).
+! Require: BND(1,J) <= BND(2,J).
+!
+! X(N) [INTENT(Out)]
+! On entry X() need not be initialized. On return,
+! X() will contain the solution N-vector.
+!
+! RNORM [INTENT(Out)]
+! The Euclidean norm of the residual vector, b - A*X.
+!
+! NSETP [INTENT(Out)]
+! Indicates the number of components of the solution
+! vector, X(), that are not at their constraint values.
+!
+! W(N) [INTENT(Out)]
+! An N-array. On return, W() will contain the dual solution
+! vector. Using Set definitions below:
+! W(J) = 0 for all j in Set P,
+! W(J) <= 0 for all j in Set Z, such that X(J) is at its
+! lower bound, and
+! W(J) >= 0 for all j in Set Z, such that X(J) is at its
+! upper bound.
+! If BND(1,J) = BND(2,J), so the variable X(J) is fixed,
+! then W(J) will have an arbitrary value.
+!
+! INDEX(N) [INTENT(Out)]
+! An INTEGER working array of size N. On exit the contents
+! of this array define the sets P, Z, and F as follows:
+!
+! INDEX(1) through INDEX(NSETP) = Set P.
+! INDEX(IZ1) through INDEX(IZ2) = Set Z.
+! INDEX(IZ2+1) through INDEX(N) = Set F.
+! IZ1 = NSETP + 1 = NPP1
+! Any of these sets may be empty. Set F is those components
+! that are constrained to a unique value by the given
+! constraints. Sets P and Z are those that are allowed a non-
+! zero range of values. Of these, set Z are those whose final
+! value is a constraint value, while set P are those whose
+! final value is not a constraint. The value of IZ2 is not returned.
+!... It is computable as the number of bounds constraining a component
+!... of X uniquely.
+!
+! IERR [INTENT(Out)]
+! Indicates status on return.
+! = 0 Solution completed.
+! = 1 M <= 0 or N <= 0
+! = 2 B(:), X(:), BND(:,:), W(:), or INDEX(:) size or shape violation.
+! = 3 Input bounds are inconsistent.
+! = 4 Exceed maximum number of iterations.
+!
+! Selected internal variables:
+! EPS [real*8]
+! Determines the relative linear dependence of a column vector
+! for a variable moved from its initial value. This is used in
+! one place with the default value EPS=EPSILON(ONE). Other
+! values, larger or smaller may be needed for some problems.
+! Library software will likely make this an optional argument.
+!
+! ITMAX [integer]
+! Set to 3*N. Maximum number of iterations permitted.
+! This is usually larger than required. Library software will
+! likely make this an optional argument.
+! ITER [integer]
+! Iteration counter.
+
+ implicit none
+ logical FIND, HITBND, FREE1, FREE2, FREE
+ integer IERR, M, N
+ integer I, IBOUND, II, IP, ITER, ITMAX, IZ, IZ1, IZ2
+ integer J, JJ, JZ, L, LBOUND, NPP1, NSETP
+
+ integer INDEX(N)
+!
+! Z() [Scratch] An M-array (automatic) of working space.
+! S() [Scratch] An N-array (automatic) of working space.
+!
+ real*8, parameter :: ZERO = 0E0, ONE=1E0, TWO = 2E0
+ real*8 :: A(M,N), B(M), S(N), X(N), W(N),
+ + Z(M), BND(2,N)
+ real*8 ALPHA, ASAVE, CC, EPS, RANGE, RNORM
+ real*8 NORM, SM, SS, T, UNORM, UP, ZTEST
+
+ CALL INITIALIZE
+!
+! The above call will set IERR.
+!
+ LOOPA: DO
+!
+! Quit on error flag, or if all coefficients are already in the
+! solution, .or. if M columns of A have been triangularized.
+!
+ IF (IERR.ne.0 .or. IZ1.gt.IZ2 .or. NSETP.ge.M) exit LOOPA
+
+ CALL SELECT_ANOTHER_COEFF_TO!_SOLVE_FOR
+!
+! See if no index was found to be moved from set Z to set P.
+! Then go to termination.
+!
+ IF ( .not. FIND ) exit LOOPA
+
+ CALL MOVE_J_FROM_SET_Z_TO_SET_P
+
+ CALL TEST_SET_P_AGAINST_CONSTRAINTS
+!
+! The above call may set IERR.
+! All coefficients in set P are strictly feasible. Loop back.
+!
+ END DO LOOPA
+
+ CALL TERMINATION
+ RETURN
+
+ CONTAINS ! These are internal subroutines.
+ SUBROUTINE INITIALIZE
+
+!*****************************************************************************80
+!
+!! INITIALIZE initializes internal data.
+!
+ IF (M.le.0 .or. N.le.0) then
+ IERR = 1
+ RETURN
+ END IF
+!
+! Check array sizes for consistency and with M and N.
+!
+ IF(SIZE(X).lt.N) THEN
+ IERR=2
+ RETURN
+ END IF
+ IF(SIZE(B).lt.M) THEN
+ IERR=2
+ RETURN
+ END IF
+ IF(SIZE(BND,1).ne.2) THEN
+ IERR=2
+ RETURN
+ END IF
+ IF(SIZE(BND,2).lt.N) THEN
+ IERR=2
+ RETURN
+ END IF
+ IF(SIZE(W).lt.N) THEN
+ IERR=2
+ RETURN
+ END IF
+ IF(SIZE(INDEX).lt.N) THEN
+ IERR=2
+ RETURN
+ END IF
+
+ IERR = 0
+!
+! The next two parameters can be changed to be optional for library software.
+!
+ EPS = EPSILON(ONE)
+ ITMAX=3*N
+ ITER=0
+!
+! Initialize the array index().
+!
+ DO I=1,N
+ INDEX(I)=I
+ END DO
+
+ IZ2=N
+ IZ1=1
+ NSETP=0
+ NPP1=1
+!
+! Begin: Loop on IZ to initialize X().
+!
+ IZ=IZ1
+ DO
+ IF (IZ.gt.IZ2 ) EXIT
+ J=INDEX(IZ)
+ IF ( BND(1,J).le.-huge(ONE)) then
+ IF (BND(2,J).ge.huge(ONE)) then
+ X(J) = ZERO
+ else
+ X(J) = min(ZERO,BND(2,J))
+ END IF
+ ELSE IF ( BND(2,J).ge.huge(ONE)) then
+ X(J) = max(ZERO,BND(1,J))
+ else
+ RANGE = BND(2,J) - BND(1,J)
+ IF ( RANGE.le.ZERO ) then
+!
+! Here X(J) is constrained to a single value.
+!
+ INDEX(IZ)=INDEX(IZ2)
+ INDEX(IZ2)=J
+ IZ=IZ-1
+ IZ2=IZ2-1
+ X(J)=BND(1,J)
+ W(J)=ZERO
+ ELSE IF ( RANGE.gt.ZERO) then
+!
+! The following statement sets X(J) to 0 if the constraint interval
+! includes 0, and otherwise sets X(J) to the endpoint of the
+! constraint interval that is closest to 0.
+!
+ X(J) = max(BND(1,J), min(BND(2,J),ZERO))
+ else
+ IERR = 3
+ RETURN
+ END IF! ( RANGE:.)
+ END IF
+
+ IF ( abs(X(J)).gt.ZERO ) then
+!
+! Change B() to reflect a nonzero starting value for X(J).
+!
+ B(1:M)=B(1:M)-A(1:M,J)*X(J)
+ END IF
+ IZ=IZ+1
+ END DO! ( IZ <= IZ2 )
+ END SUBROUTINE! ( INITIALIZE )
+
+ SUBROUTINE SELECT_ANOTHER_COEFF_TO!_SOLVE_FOR
+
+!*****************************************************************************80
+!
+!! SELECT_ANOTHER_COEFF_TO!_SOLVE_FOR
+!
+! 1. Search through set z for a new coefficient to solve for.
+! First select a candidate that is either an unconstrained
+! coefficient or else a constrained coefficient that has room
+! to move in the direction consistent with the sign of its dual
+! vector component. Components of the dual (negative gradient)
+! vector will be computed as needed.
+! 2. For each candidate start the transformation to bring this
+! candidate into the triangle, and then do two tests: Test size
+! of new diagonal value to avoid extreme ill-conditioning, and
+! the value of this new coefficient to be sure it moved in the
+! expected direction.
+! 3. If some coefficient passes all these conditions, set FIND = true,
+! The index of the selected coefficient is J = INDEX(IZ).
+! 4. If no coefficient is selected, set FIND = false.
+!
+ FIND = .FALSE.
+ DO IZ=IZ1,IZ2
+ J=INDEX(IZ)
+!
+! Set FREE1 = true if X(J) is not at the left end-point of its
+! constraint region.
+! Set FREE2 = true if X(J) is not at the right end-point of its
+! constraint region.
+! Set FREE = true if X(J) is not at either end-point of its
+! constraint region.
+!
+ FREE1 = X(J).gt.BND(1,J)
+ FREE2 = X(J).lt.BND(2,J)
+ FREE = FREE1 .and. FREE2
+
+ IF ( FREE ) then
+ CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+ else
+!
+! Compute dual coefficient W(J).
+!
+ W(J)=dot_product(A(NPP1:M,J),B(NPP1:M))

Use BLAS *DOT instead of intrinsic fortran 90 dot_product

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+!
+! Can X(J) move in the direction indicated by the sign of W(J)?
+!
+ IF ( W(J).lt.ZERO ) then
+ IF ( FREE1 )
+ + CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+ ELSE IF ( W(J).gt.ZERO ) then
+ IF ( FREE2 )
+ + CALL TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+ END IF
+ END IF
+ IF ( FIND ) RETURN
+ END DO! IZ
+ END SUBROUTINE! ( SELECT ANOTHER COEF TO SOLVE FOR )
+
+ SUBROUTINE TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+
+!*****************************************************************************80
+!
+!! TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+!
+! The sign of W(J) is OK for J to be moved to set P.
+! Begin the transformation and check new diagonal element to avoid
+! near linear dependence.
+!
+ ASAVE=A(NPP1,J)
+ call HTC (NPP1, A(1:M,J), UP)
+ UNORM = NRM2(A(1:NSETP,J))
+ IF ( abs(A(NPP1,J)).gt.EPS * UNORM) then
+!
+! Column J is sufficiently independent. Copy b into Z, update Z.
+!
+ Z(1:M)=B(1:M)
+!
+! Compute product of transormation and updated right-hand side.
+!
+ NORM=A(NPP1,J); A(NPP1,J)=UP
+ IF(ABS(NORM).gt.ZERO) THEN
+ SM=DOT_PRODUCT(A(NPP1:M,J)/NORM, Z(NPP1:M))/UP
+ Z(NPP1:M)=Z(NPP1:M)+SM*A(NPP1:M,J)
+ A(NPP1,J)=NORM
+ END IF
+
+ IF (abs(X(J)).gt.ZERO) Z(1:NPP1)=Z(1:NPP1)+A(1:NPP1,J)*X(J)
+!
+! Adjust Z() as though X(J) had been reset to zero.
+!
+ IF ( FREE ) then
+ FIND = .TRUE.
+ else
+!
+! Solve for ZTEST ( proposed new value for X(J) ).
+! Then set FIND to indicate if ZTEST has moved away from X(J) in
+! the expected direction indicated by the sign of W(J).
+!
+ ZTEST=Z(NPP1)/A(NPP1,J)
+ FIND = ( W(J).lt.ZERO .and. ZTEST.lt.X(J) ) .or.
+ + ( W(J).gt.ZERO .and. ZTEST.gt.X(J) )
+ END IF
+ END IF
+!
+! If J was not accepted to be moved from set Z to set P,
+! restore A(NNP1,J). Failing these tests may mean the computed
+! sign of W(J) is suspect, so here we set W(J) = 0. This will
+! not affect subsequent computation, but cleans up the W() array.
+!
+ IF ( .not. FIND ) then
+ A(NPP1,J)=ASAVE
+ W(J)=ZERO
+ END IF! ( .not. FIND )
+ END SUBROUTINE !TEST_COEF_J_FOR_DIAG!_ELT_AND_DIRECTION_OF_CHANGE
+
+ SUBROUTINE MOVE_J_FROM_SET_Z_TO_SET_P
+
+!*****************************************************************************80
+!
+!! MOVE_J_FROM_SET_Z_TO_SET_P
+!
+! The index J=index(IZ) has been selected to be moved from
+! set Z to set P. Z() contains the old B() adjusted as though X(J) = 0.
+! A(*,J) contains the new Householder transformation vector.
+!
+ B(1:M)=Z(1:M)
+
+ INDEX(IZ)=INDEX(IZ1)
+ INDEX(IZ1)=J
+ IZ1=IZ1+1
+ NSETP=NPP1
+ NPP1=NPP1+1
+!
+! The following loop can be null or not required.
+!
+ NORM=A(NSETP,J); A(NSETP,J)=UP
+ IF(ABS(NORM).gt.ZERO) THEN
+ DO JZ=IZ1,IZ2
+ JJ=INDEX(JZ)
+ SM=DOT_PRODUCT(A(NSETP:M,J)/NORM, A(NSETP:M,JJ))/UP
+ A(NSETP:M,JJ)=A(NSETP:M,JJ)+SM*A(NSETP:M,J)
+ END DO
+ A(NSETP,J)=NORM
+ END IF
+!
+! The following loop can be null.
+!
+ DO L=NPP1,M
+ A(L,J)=ZERO
+ END DO! L
+
+ W(J)=ZERO
+!
+! Solve the triangular system. Store this solution temporarily in Z().
+!
+ DO I = NSETP, 1, -1
+ IF (I.ne.NSETP) Z(1:I)=Z(1:I)-A(1:I,II)*Z(I+1)
+ II=INDEX(I)
+ Z(I)=Z(I)/A(I,II)
+ END DO
+ END SUBROUTINE! ( MOVE J FROM SET Z TO SET P )
+
+ SUBROUTINE TEST_SET_P_AGAINST_CONSTRAINTS
+
+!*****************************************************************************80
+!
+!! TEST_SET_P_AGAINST_CONSTRAINTS
+!
+ LOOPB: DO
+!
+! The solution obtained by solving the current set P is in the array Z().
+!
+ ITER=ITER+1
+ IF (ITER.gt.ITMAX) then
+ IERR = 4
+ exit LOOPB
+ END IF
+
+ CALL SEE_IF_ALL_CONSTRAINED_COEFFS!_ARE_FEASIBLE
+!
+! The above call sets HITBND. If HITBND = true then it also sets
+! ALPHA, JJ, and IBOUND.
+!
+ IF ( .not. HITBND ) exit LOOPB
+!
+! Here ALPHA will be between 0 and 1 for interpolation
+! between the old X() and the new Z().
+!
+ DO IP=1,NSETP
+ L=INDEX(IP)
+ X(L)=X(L)+ALPHA*(Z(IP)-X(L))
+ END DO
+
+ I=INDEX(JJ)
+!
+! Note: The exit test is done at the end of the loop, so the loop
+! will always be executed at least once.
+!
+ DO
+!
+! Modify A(*,*), B(*) and the index arrays to move coefficient I
+! from set P to set Z.
+!
+ CALL MOVE_COEF_I_FROM_SET_P_TO_SET_Z
+
+ IF (NSETP.le.0) exit LOOPB
+!
+! See if the remaining coefficients in set P are feasible. They should
+! be because of the way ALPHA was determined. If any are infeasible
+! it is due to round-off error. Any that are infeasible or on a boundary
+! will be set to the boundary value and moved from set P to set Z.
+!
+ IBOUND = 0
+ DO JJ=1,NSETP
+ I=INDEX(JJ)
+ IF ( X(I).le.BND(1,I)) then
+ IBOUND=1
+ EXIT
+ ELSE IF ( X(I).ge.BND(2,I)) then
+ IBOUND=2
+ EXIT
+ END IF
+ END DO
+ IF (IBOUND.le.0) EXIT
+ END DO
+!
+! Copy B( ) into Z( ). Then solve again and loop back.
+!
+ Z(1:M)=B(1:M)
+
+ DO I = NSETP, 1, -1
+ IF (I.ne.NSETP) Z(1:I)=Z(1:I)-A(1:I,II)*Z(I+1)

Vectorized expressions are invalid in Fortran 77.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ II=INDEX(I)
+ Z(I)=Z(I)/A(I,II)
+ END DO
+ END DO LOOPB
+!
+! The following loop can be null.
+!
+ DO IP=1,NSETP
+ I=INDEX(IP)
+ X(I)=Z(IP)
+ END DO
+
+ END SUBROUTINE! ( TEST SET P AGAINST CONSTRAINTS)
+
+ SUBROUTINE SEE_IF_ALL_CONSTRAINED_COEFFS!_ARE_FEASIBLE
+
+!*****************************************************************************80
+!
+!! SEE_IF_ALL_CONSTRAINED_COEFFS!_ARE_FEASIBLE
+!
+! See if each coefficient in set P is strictly interior to its constraint
+! region.
+! If so, set HITBND = false.
+! If not, set HITBND = true, and also set ALPHA, JJ, and IBOUND.
+! Then ALPHA will satisfy 0. < ALPHA <= 1.
+!
+ ALPHA=TWO
+ DO IP=1,NSETP
+ L=INDEX(IP)
+ IF (Z(IP).le.BND(1,L)) then
+!
+! Z(IP) HITS LOWER BOUND
+!
+ LBOUND=1
+ ELSE IF (Z(IP).ge.BND(2,L)) then
+!
+! Z(IP) HITS UPPER BOUND
+!
+ LBOUND=2
+ else
+ LBOUND = 0
+ END IF
+
+ IF ( LBOUND.ne.0 ) then

Using name of intrinsic function as variable name.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ T=(BND(LBOUND,L)-X(L))/(Z(IP)-X(L))
+ IF (ALPHA.gt.T) then
+ ALPHA=T
+ JJ=IP
+ IBOUND=LBOUND
+ END IF! ( LBOUND )
+ END IF! ( ALPHA > T )
+ END DO
+ HITBND = abs(ALPHA - TWO).gt.ZERO
+ END SUBROUTINE!( SEE IF ALL CONSTRAINED COEFFS ARE FEASIBLE )
+
+ SUBROUTINE MOVE_COEF_I_FROM_SET_P_TO_SET_Z
+
+!*****************************************************************************80
+!
+!! MOVE_COEF_I_FROM_SET_P_TO_SET_Z
+!
+ X(I)=BND(IBOUND,I)
+ IF (abs(X(I)).gt.ZERO .and. JJ.gt.0) B(1:JJ)=B(1:JJ)-A(1:JJ,I)*X(I)
+!
+! The following loop can be null.
+!
+ DO J = JJ+1, NSETP
+ II=INDEX(J)
+ INDEX(J-1)=II
+ call ROTG (A(J-1,II),A(J,II),CC,SS)
+ SM=A(J-1,II)
+!
+! The plane rotation is applied to two rows of A and the right-hand
+! side. One row is moved to the scratch array S and then the updates
+! are computed. The intent is for array operations to be performed
+! and minimal extra data movement. One extra rotation is applied
+! to column II in this approach.
+!
+ S=A(J-1,1:N); A(J-1,1:N)=CC*S+SS*A(J,1:N); A(J,1:N)=
+ + CC*A(J,1:N)-SS*S
+ A(J-1,II)=SM; A(J,II)=ZERO
+ SM=B(J-1); B(J-1)=CC*SM+SS*B(J); B(J)=CC*B(J)-SS*SM
+ END DO
+
+ NPP1=NSETP
+ NSETP=NSETP-1
+ IZ1=IZ1-1
+ INDEX(IZ1)=I
+ END SUBROUTINE! ( MOVE COEF I FROM SET P TO SET Z )
+
+ SUBROUTINE TERMINATION
+
+!*****************************************************************************80
+!
+!! TERMINATION
+!
+ IF (IERR.le.0) then
+!
+! Compute the norm of the residual vector.
+!
+ SM=ZERO
+ IF (NPP1.le.M) then
+ SM=NRM2(B(NPP1:M))
+ else
+ W(1:N)=ZERO
+ END IF
+ RNORM=SM
+ END IF! ( IERR...)
+ END SUBROUTINE! ( TERMINATION )
+
+ SUBROUTINE ROTG(SA,SB,C,S)
+
+!*****************************************************************************80
+!
+!! ROTG performs a Givens rotation.
+!
+ REAL*8 SA,SB,C,S,ROE,SCALE,R
+
+ ROE = SB
+ IF( ABS(SA) .GT. ABS(SB) ) ROE = SA
+ SCALE = ABS(SA) + ABS(SB)
+ IF( SCALE .LE. ZERO ) THEN
+ C = ONE
+ S = ZERO
+ RETURN
+ END IF
+ R = SCALE*SQRT((SA/SCALE)**2 + (SB/SCALE)**2)
+ IF(ROE.lt.ZERO)R=-R
+ C = SA/R
+ S = SB/R
+ SA = R
+ RETURN
+ END SUBROUTINE !ROTG
+
+ REAL*8 FUNCTION NRM2 (X)
+
+!*****************************************************************************80
+!
+!! NRM2 returns the Euclidean norm of a vector.
+!
+! NRM2 := sqrt( x'*x )
+!
+ REAL*8 ABSXI, X(:), NORM, SCALE, SSQ
+ INTEGER N, IX
+ N=SIZE(X)
+ IF( N.lt.1)THEN
+ NORM = ZERO
+ ELSE IF( N.eq.1 )THEN
+ NORM = ABS( X( 1 ) )
+ ELSE
+ SCALE = ZERO
+ SSQ = ONE
+
+ DO IX = 1, N
+ ABSXI = ABS( X( IX ) )
+ IF(ABSXI.gt.ZERO )THEN
+ IF( SCALE.lt.ABSXI )THEN
+ SSQ = ONE + SSQ*( SCALE/ABSXI )**2
+ SCALE = ABSXI
+ ELSE
+ SSQ = SSQ + ( ABSXI/SCALE )**2
+ END IF
+ END IF
+ END DO
+ NORM = SCALE * SQRT( SSQ )
+ END IF
+
+ NRM2 = NORM
+ RETURN
+ END FUNCTION
+
+ SUBROUTINE HTC (P, U, UP)
+
+!*****************************************************************************80
+!
+!! HTC constructs a Householder transformation.
+!
+ INTEGER P
+ REAL*8 U(:)
+ REAL*8 UP, VNORM
+ VNORM=NRM2(U(P:SIZE(U)))
+ IF(U(P).gt.ZERO) VNORM=-VNORM
+ UP=U(P)-VNORM
+ U(P)=VNORM
+ END SUBROUTINE ! HTC
+
+ END SUBROUTINE ! BVLS
View
23 scipy/optimize/bvls/bvls.pyf
@@ -0,0 +1,23 @@
+! -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+python module _bvlslib ! in
+ interface ! in :_bvlslib
+ subroutine bvls(a,b,m,n,bnd,x,rnorm,nsetp,w,index_bn,ierr) ! in :_bvlslib:bvls.f
+ real*8 dimension(m,n),intent(inout) :: a
+ real*8 dimension(m),depend(m),intent(inout) :: b
+ integer, optional,check(shape(a,0)==m),depend(a),intent(hide) :: m=shape(a,0)
+ integer, optional,check(shape(a,1)==n),depend(a),intent(hide) :: n=shape(a,1)
+ real*8 dimension(2,n),depend(n),intent(in) :: bnd
+ real*8 dimension(n),depend(n),intent(inout) :: x
+ real*8,intent(out) :: rnorm
+ integer,intent(out) :: nsetp
+ real*8 dimension(n),depend(n),intent(in) :: w
+ integer dimension(n),depend(n),intent(in) :: index_bn
+ integer,intent(out) :: ierr
+ end subroutine bvls
+ end interface
+end python module _bvls
+
+! This file was auto-generated with f2py (version:2).
+! See http://cens.ioc.ee/projects/f2py2e/
View
4 scipy/optimize/setup.py
@@ -53,6 +53,10 @@ def configuration(parent_package='',top_path=None):
config.add_extension('_nnls', sources=[join('nnls', x)
for x in ["nnls.f","nnls.pyf"]])
+ sources = ["bvls.f", "bvls.pyf"]
+ config.add_extension('_bvlslib', sources=[join('bvls', x)
+ for x in sources])
+
config.add_data_dir('tests')
config.add_data_dir('benchmarks')
return config
View
30 scipy/optimize/tests/test_bvls.py
@@ -0,0 +1,30 @@
+from scipy.optimize import bounded_lstsq
+import numpy as np
+from numpy.linalg import norm, lstsq
+from numpy.testing import (TestCase, assert_array_almost_equal, assert_,
+ assert_allclose)
+
+
+class TestBVLS(TestCase):
+ def test_bvls(self):
+ np.random.seed(1)
+ X = np.random.randn(25, 5)
+ b = np.arange(1, 6.)
+ y = np.dot(X, b) + np.random.randn(25)
+ bb, res, nsetp = bounded_lstsq(X, y, bounds=[(0, None), (1, None),
+ (2, None), (3, None),
+ (4, None)])
+ expected, _, _, _ = lstsq(X, y)
+ assert_array_almost_equal(bb, expected)
+ assert_allclose(norm(np.dot(X, bb)-y), res)
+
+ bb, res, nsetp = bounded_lstsq(X, y, bounds=[(None, None),
+ (None, None),
+ (None, None),
+ (1, 2),
+ (2, 3)])
+ assert_(bb[3] == 2)
+ assert_(bb[4] == 3)
+
+ # smoke test
+ bb, res, nsetp = bounded_lstsq(X, y, bounds=None)
Something went wrong with that request. Please try again.