Mailman 3 July 2010 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck 24 Nov '13

24 Nov '13

I'm currently implementing the Sequential Least Squares Quadratic Programming (SLSQP) optimizer, by Dieter Kraft, for use in Scipy. The Fortran code being wrapped with F2PY is here: http://www.netlib.org/toms/733 (its use within Scipy has been cleared) SLSQP provides for bounds on the independent variables, as well as equality and inequality constraint functions, which is a capability that doesn't exist in scipy.optimize. Currently the code works, although the constraint normals are being generated by approximation. I'm working on a way to pass these in. I think the most elegant way will be a single function that returns the matrix of constraint normals. For a demonstration of what the code can do, here is an optimization of f(x,y) = 2xy + 2x - x**2 - 2y**2 Example 14.2 in Chapra & Canale gives the maximum as x=2.0, y=1.0. The unbounded optimization tests find this solution. As expected, its faster when derivatives are provided rather than approximated. Unbounded optimization. Derivatives approximated. Elapsed time: 1.45792961121 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] Unbounded optimization. Derivatives provided. Elapsed time: 1.03211402893 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] The following example uses an equality constraint to find the optimal when x=y. Bound optimization. Derivatives approximated. Elapsed time: 1.384973526 ms Results [[0.99999996845920858, 0.99999996845920858], -0.99999999999999889, 4, 0, 'Optimization terminated successfully.'] I've tried to conform to the syntax used by the other optimizers in scipy.optimize. The function definition and doc string are below. If anyone is interested in testing it out, let me know. def fmin_slsqp( func, x0 , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None, bounds = [], fprime = None, fprime_cons=None,args = (), iter = 100, acc = 1.0E-6, iprint = 1, full_output = 0, epsilon = _epsilon ): """ Minimize a function using Sequential Least SQuares Programming Description: Python interface function for the SLSQP Optimization subroutine originally implemented by Dieter Kraft. Inputs: func - Objective function (in the form func(x, *args)) x0 - Initial guess for the independent variable(s). eqcons - A list of functions of length n such that eqcons[j](x0,*args) == 0.0 in a successfully optimized problem f_eqcons - A function of the form f_eqcons(x, *args) that returns an array in which each element must equal 0.0 in a successfully optimized problem. If f_eqcons is specified, eqcons is ignored. ieqcons - A list of functions of length n such that ieqcons[j](x0,*args) >= 0.0 in a successfully optimized problem f_ieqcons - A function of the form f_ieqcons(x0, *args) that returns an array in which each element must be greater or equal to 0.0 in a successfully optimized problem. If f_ieqcons is specified, ieqcons is ignored. bounds - A list of tuples specifying the lower and upper bound for each independent variable [ (xl0, xu0), (xl1, xu1), ...] fprime - A function that evaluates the partial derivatives of func fprime_cons - A function of the form f(x, *args) that returns the m by n array of constraint normals. If not provided, the normals will be approximated. Equality constraint normals precede inequality constraint normals. args - A sequence of additional arguments passed to func and fprime iter - The maximum number of iterations (int) acc - Requested accuracy (float) iprint - The verbosity of fmin_slsqp. iprint <= 0 : Silent operation iprint == 1 : Print summary upon completion (default) iprint >= 2 : Print status of each iterate and summary full_output - If 0, return only the minimizer of func (default). Otherwise, output final objective function and summary information. epsilon - The step size for finite-difference derivative estimates. Outputs: ( x, { fx, gnorm, its, imode, smode }) x - The final minimizer of func. fx - The final value of the objective function. its - The number of iterations. imode - The exit mode from the optimizer, as an integer. smode - A string describing the exit mode from the optimizer. Exit modes are defined as follows: -1 : Gradient evaluation required (g & a) 0 : Optimization terminated successfully. 1 : Function evaluation required (f & c) 2 : Number of equality constraints larger than number of independent variables 3 : More than 3*n iterations in LSQ subproblem 4 : Inequality constraints incompatible 5 : Singular matrix E in LSQ subproblem 6 : Singular matrix C in LSQ subproblem 7 : Rank-deficient equality constraint subproblem HFTI 8 : Positive directional derivative for linesearch 9 : Iteration limit exceeded -- - Rob Falck

4 16

scipy error compiling csr_wrap under Python 2.7
by Bruce Southey 27 Sep '10

27 Sep '10

Hi, I failed to get the scipy 0.8 rc1 and SVN to build under Python2.7 with numpy '2.0.0.dev8469' and gcc version 4.4.4 20100503 (Red Hat 4.4.4-2) (GCC) . But scipy does compile with Python 2.6. building 'scipy.sparse.sparsetools._csr' extension compiling C++ sources C compiler: g++ -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -fPIC compile options: '-I/usr/local/lib/python2.7/site-packages/numpy/core/include -I/usr/local/include/python2.7 -c' g++: scipy/sparse/sparsetools/csr_wrap.cxx scipy/sparse/sparsetools/csr_wrap.cxx: In function ‘int require_size(PyArrayObject*, npy_intp*, int)’: scipy/sparse/sparsetools/csr_wrap.cxx:2910: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx:2917: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx: In function ‘int require_size(PyArrayObject*, npy_intp*, int)’: scipy/sparse/sparsetools/csr_wrap.cxx:2910: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2910: warning: too many arguments for format scipy/sparse/sparsetools/csr_wrap.cxx:2917: error: expected ‘)’ before ‘PRIdPTR’ scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: spurious trailing ‘%’ in format scipy/sparse/sparsetools/csr_wrap.cxx:2917: warning: too many arguments for format error: Command "g++ -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -fPIC -I/usr/local/lib/python2.7/site-packages/numpy/core/include -I/usr/local/include/python2.7 -c scipy/sparse/sparsetools/csr_wrap.cxx -o build/temp.linux-x86_64-2.7/scipy/sparse/sparsetools/csr_wrap.o" failed with exit status 1 The error is possibly related the variable 'NPY_INTP_FMT' (defined in numpy/core/include/numpy/ndarraytypes.h). A grep shows these files: scipy/sparse/sparsetools/coo_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/coo_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/bsr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/bsr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/numpy.i: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/numpy.i: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/csr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/csr_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/csc_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/csc_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); scipy/sparse/sparsetools/dia_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", size[i]); scipy/sparse/sparsetools/dia_wrap.cxx: sprintf(s,"%" NPY_INTP_FMT ",", array_size(ary,i)); Note that the scipy/sparse/sparsetools/numpy.i contains tabs in some places instead of spaces. The ndarraytypes.h is present in /usr/local/lib/python2.7/site-packages/numpy/core/include/numpy/ Bruce

5 10

N-dimensional interpolation
by Pauli Virtanen 06 Sep '10

06 Sep '10

Hi all, I took the Qhull by the horns, and wrote a straightforward `griddata` implementation for working in N-D: http://github.com/pv/scipy-work/tree/qhull http://github.com/pv/scipy-work/blob/qhull/scipy/interpolate/qhull.pyx Sample (4-D): ------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import griddata def f(x, y, z, u): return np.cos(x + np.cos(y + np.cos(z + np.cos(u)))) points = np.random.rand(2000, 4) values = f(*points.T) p = np.linspace(0, 1, 1500) z = griddata(points, values, np.c_[p, p, p, p], method='linear') plt.plot(p, z, '-', p, f(p, p, p, p), '-') plt.show() ------------------------------------------------------------- It performs the N-D Delaunay tesselation with Qhull, and after that does linear barycentric interpolation on the simplex containing each point. (The nearest-neighbor "interpolation" is, on the other hand, implemented using scipy.spatial.KDTree.) I ended up writing some custom simplex-walking code for locating the simplex containing the points -- this is much faster than a brute-force search. However, it is slightly different from what `qh_findbestfacet` does. [If you're a computational geometry expert, see below...] Speed-wise, Qhull appears to lose in 2-D to `scikits.delaunay` by a factor of about 10 in triangulation in my tests; this, however, includes the time taken by computing the barycentric coordinate transforms. Interpolation, however, seems to be faster. *** I think this would be a nice addition to `scipy.optimize`. I'd like to get it in for 0.9. Later on, we can specialize the `griddata` function to perform better in low dimensions (ndim=2, 3) and add more interpolation methods there; for example natural neighbors, interpolating splines, etc. In ndim > 3, the `griddata` is however already feature-equivalent to Octave and MATLAB. Comments? Pauli PS. A question for computational geometry experts: `qh_findbestfacet` finds the facet on the lower convex hull whose oriented hyperplane distance to the point lifted on the paraboloid is maximal --- however, this does not always give the simplex containing the point in the projected Delaunay tesselation. There's a counterexample in qhull.pyx:_find_simplex docstring. On the other hand, Qhull documentation claims that this is the way to locate the Delaunay simplex containing a point. What gives? It's clear, however, that if a simplex contains the point, then the hyperplane distance is positive. So currently, I just walk around positive-distance simplices checking the inside condition for each and hoping for the best. If nothing is found, I just fall back to brute force search. This seems to work well in practice. However, if a point is outside the triangulation (but inside the bounding box), I have to fall back to a brute-force search. Is there a better option here? *** Another sample (2-D, object-oriented interface): ------------------------------------------------------------- import time import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import LinearNDInterpolator # some random data x = np.array([(0,0), (0,1), (1, 1), (1, 0)], dtype=np.double) np.random.seed(0) xr = np.random.rand(200*200, 2).astype(np.double) x = np.r_[xr, x] y = np.arange(x.shape[0], dtype=np.double) # evaluate on a grid nn = 500 xx = np.linspace(-0.1, 1.1, nn) yy = np.linspace(-0.1, 1.1, nn) xx, yy = np.broadcast_arrays(xx[:,None], yy[None,:]) xi = np.array([xx, yy]).transpose(1,2,0).copy() start = time.time() ip = LinearNDInterpolator(x, y) print "Triangulation", time.time() - start start = time.time() zi = ip(xi) print "Interpolation", time.time() - start from scikits.delaunay import Triangulation start = time.time() tri2 = Triangulation(x[:,0], x[:,1]) print "scikits.delaunay triangulation", time.time() - start start = time.time() ip2 = tri2.linear_interpolator(y) zi2 = ip2[-0.1:1.1:500j,-0.1:1.1:500j].T print "scikits.delaunay interpolation", time.time() - start print "rel-difference", np.nanmax(abs(zi - zi2))/np.nanmax(np.abs(zi)) plt.figure() plt.imshow(zi) plt.clim(y.min(), y.max()) plt.colorbar() plt.figure() plt.imshow(zi2) plt.clim(y.min(), y.max()) plt.colorbar() plt.show() ------------------------------------------------------------- -- Pauli Virtanen

5 8

Getting Scipy's weave to work reliably on Windows
by Chris Ball 24 Aug '10

24 Aug '10

Hi, While testing Scipy's weave on several different Windows installations, I came across some problems with spaces in paths that often prevent weave from working. I can see a change that could probably get weave working on most Windows installations, but it is a quick hack. Someone knowledgeable about distutils (and numpy.distutils?) might be able to help me fix this properly. Below I describe three common problems with weave on Windows, in the hope that this information helps others, or allows someone to suggest how to fix the spaces-in- paths problem properly. I think there are three common problems that stop weave from working on Windows. The first is not having a C compiler. Both Python(x,y) and EPD provide a C compiler that seems to work fine, which is great! The second problem is that if weave is installed to a location with a space in the path, linking fails. There is already a scipy bug report about this (http://projects.scipy.org/scipy/ticket/809). I've just commented on that report, saying the problem appears to be with distutils, and there is already a Python bug report about it (http://bugs.python.org/issue4508). Maybe someone could close this scipy bug, or link it to the Python one somehow? In any case, when using Python(x,y) or EPD, this bug will not show up if the default installation locations are accepted. So, that's also good news! The third problem is that if the Windows user name has a space in it (which in my experience is quite common), compilation fails. Weave uses the user name to create a path for its "intermediate" and "compiled" files. When the compilation command is issued, the path with the space in it is also not quoted. Presumably that is another error in distutils (or numpy.distutils)? Unfortunately I wasn't able to pinpoint what function is failing to quote strings properly, because I couldn't figure out the chain that leads to the compiler being called. However, I can avoid the problem by removing spaces from the user name in weave itself (catalog.py): def whoami(): """return a string identifying the user.""" return (os.environ.get("USER") or os.environ.get("USERNAME") or "unknown").replace(" ","") (where I have added .replace(" ","") to the existing code). I realize this isn't the right solution, so if someone could help to guide me to the point where quoting should occur, that would be very helpful. Otherwise, is there any chance of applying a hack like this so weave can work reliably on Windows? Thanks, Chris

2 1

Would it be acceptable to give up subpackage build ?
by David Cournapeau 31 Jul '10

31 Jul '10

Hi there, I have made enough progress on bento recently to start tackling the meat of the issue: building numpy and scipy with it. Numpy goes nicely, but scipy causes issues because of the following constraing: in scipy, each package with a setup.py is a self-contained package, i.e. you can just do python setup.py build in scipy/lib/blas, and it will run. Is it acceptable to give up on this in the long run ? I don't see much advantages in that feature, and giving it up has several nice advantages: - one can have a global configure stage (i.e. check blas/lapack only once instead of many times as done currently, although the way numpy.distutils has near 0 cost since it only stat files) - I would expect the build to be faster and more accurate - it makes my life easier for bento Note that this has no incidence on the existing distutils-based infrastructure, cheers, David

3 2

Proposed enhancement: pure Python implementation of LP simplex method (two-phase)
by Enzo Michelangeli 30 Jul '10

30 Jul '10

I just submitted the following enhancement ticket: http://projects.scipy.org/scipy/ticket/1252 Enzo

4 6

I'm getting a "OperationalError: database is locked"
by David Goldsmith 29 Jul '10

29 Jul '10

when trying to file a new scipy ticket; I am logged in. DG

2 1

Good-bye, sort of
by David Goldsmith 29 Jul '10

29 Jul '10

July 31 is my last day as paid NumPy/SciPy docs editor/marathon coordinator; I'm moving on to a programming data analyst position for the Marine Waters Data Management System Project of the Environmental Assessment Program in the Washington State Department of Ecology; unfortunately, they need me to program in MATLAB, but I plan on continuing to use NumPy, etc. for my fractals, etc. I also intend to continue to contribute to this community as a volunteer as my time permits, but you probably will be seeing less of me on the lists. Also, I won't be hosting a Marathon Skypecon this Friday (but of course others are free to organize one if they want to.) Thanks for all your help and support, David Goldsmith Olympia, WA

8 7

numpy svn broken NPY_COPY_PYOBJECT_PTR
by Nils Wagner 29 Jul '10

29 Jul '10

Hi all, I have installed the most recent svn version of numpy. However if I try to import numpy it fails with >>> import numpy Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/__init__.py", line 136, in <module> import add_newdocs File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/add_newdocs.py", line 9, in <module> from numpy.lib import add_newdoc File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/lib/__init__.py", line 4, in <module> from type_check import * File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/lib/type_check.py", line 8, in <module> import numpy.core.numeric as _nx File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/core/__init__.py", line 5, in <module> import multiarray ImportError: /data/home/nwagner/local/lib/python2.5/site-packages/numpy/core/multiarray.so: undefined symbol: NPY_COPY_PYOBJECT_PTR I cannot install the most recent svn version of scipy python setup.py install --prefix=$HOME/local Traceback (most recent call last): File "setup.py", line 160, in <module> setup_package() File "setup.py", line 127, in setup_package from numpy.distutils.core import setup File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/__init__.py", line 136, in <module> import add_newdocs File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/add_newdocs.py", line 9, in <module> from numpy.lib import add_newdoc File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/lib/__init__.py", line 4, in <module> from type_check import * File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/lib/type_check.py", line 8, in <module> import numpy.core.numeric as _nx File "/data/home/nwagner/local/lib/python2.5/site-packages/numpy/core/__init__.py", line 5, in <module> import multiarray ImportError: /data/home/nwagner/local/lib/python2.5/site-packages/numpy/core/multiarray.so: undefined symbol: NPY_COPY_PYOBJECT_PTR Any idea ? Nils

3 3

severe bug in scipy.integrate.quadrature?
by Thomas Robitaille 28 Jul '10

28 Jul '10

Hi, I think there is a mistake in scipy.integrate.quadrature that can result in false convergence of the integral. The code in question is the following: err = 100. val = err n = 1 vfunc = vectorize1(func, args, vec_func=vec_func) while (err > tol) and (n < maxiter): newval = fixed_quad(vfunc, a, b, (), n)[0] err = abs(newval-val) val = newval n = n + 1 The error inside the loop is the absolute error, so it looks like the tolerance to specify is the absolute tolerance, not relative. However, the original err is set to 100, so if I compute an integral that has a result of 1e20, I might want a tolerance of say 1e15, in which case the loop will be skipped, because err is not greater than tol. In its present form, I think that err should be initially set to infinity. In addition, I think that the tolerance to specify should be the relative one, not the absolute, because most of the time, I know I want my integrals calculated to a precision of say 1e-5, but I don't know in advance the answer! I think a relative tolerance is what users expect at any rate. Can any experts comment on this? Thanks, Thomas

2 1