Mailman 3 October 2009 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck Nov. 24, 2013

Nov. 24, 2013

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

Ball Tree class
by Jake VanderPlas Dec. 3, 2009

Dec. 3, 2009

Hello, I've been using scipy.spatial.KDTree for my research, and found it very useful. Recently, though, my datasets have been getting larger - upwards of 10,000 points in 1000-2000 dimensions. The KDTree starts getting slow with such a large dimensionality. I've addressed this by writing a C++ Ball Tree code, and wrapping it using swig and numpy.i. I've been wanting to begin contributing to the scipy project, and I think this would be a great place to start. I'd like to begin the process of adding this to the scipy.spatial package. A few questions: Is it preferable to have an implementation in C rather than C++? Cython, swig, or hand-wrapped code? cKDTree is written in Cython, with C. Should I stick to that convention to maintain uniformity in the scipy.spatial package? Please let me know what you think -Jake

4 3

ndimage - docfiller and output_type
by Ralf Gommers Nov. 16, 2009

Nov. 16, 2009

Hi all, While working on the ndimage docs, I saw that the docfiller decorator in ndimage.filters does not work well with the doc wiki. Basically, in the doc wiki you see the docstring with the decorator applied. Then if you make an edit, all the %(param) markup has been replaced. Compare the `convolve` and `convolve1d` docstrings in trunk to see what happens. Some options to fix this: 1. remove the filters docstrings from the doc wiki 2. make the doc wiki understand the decorator 3. remove the decorator, and do some copy-paste'ing of parameter descriptions 4. remove the decorator, and refer to a single description of those parameters outside the docstring. I would have a preference for 3 or 4, mostly because it's easier and more robust. Would you like a patch for either of those options? The second thing I came across is the output_type keyword in ndimage.interpolation. Specifying anything other than None for that keyword results in a RuntimeError. The code responsible for that is in ndimage._ni_support.py : def _get_output(output, input, output_type = None, shape = None): if output_type is not None: msg = "'output_type' argument is deprecated." msg += " Assign type to 'output' instead." raise RuntimeError, msg <6 more lines that never execute> Is there any reason to keep that keyword around? I can write a patch for that as well. Cheers, Ralf

6 32

Continuous Wavelet Transform for SciPy
by Sean Arms Nov. 12, 2009

Nov. 12, 2009

Greetings! My name is Sean Arms and I'm a graduate student at the University of Oklahoma in the School of Meteorology. As part of my PhD research, I'm studying coherent structures in atmospheric boundary layer turbulence, primarily using in-situ observations and, secondarily, Direct Numerical Simulation (DNS) output. One common approach for detecting coherent structures in observational datasets relies on the use of the global wavelet power spectrum as estimated from a continuous wavelet transform (CWT). I know SciPy has a DWT impementation, and I've already been in contact with Filip. He recommeded that I post my code in hopes that it would add some momentum to the python-cwt development and create some feedback (I'm currently looking for a good place to post my code). I've implemented the CWT using pure python (that is, I did not write any C extension code myself - nothing to build), along with one mother wavelet (second derivitive of a Gaussian, or the Mexican Hat) - I'll be adding more Mother wavelets as I go along. I've made it a point to (try to) design my MotherWavelet class to be easily extendable. I'm working on documentation and a few tests at the moment, but so far my code compares well with other wavelet routines. The point of this email is to introduce myself and let the SciPy dev community know that I am willing to help develop CWT support for SciPy - I'll already be doing the work for my research, so I might as well put in the extra effort to make is usable by the larger community! Cheers! Sean Arms Graduate Research Assistant School of Meteorology University of Oklahoma

4 7

Lambert W implementation
by Yosef Meller Nov. 7, 2009

Nov. 7, 2009

Hello, I searched through scipy.special and I saw no implementation of the Lambert W function. I'd like to add one, but before I do I thought I'd ask: 1) Did I just miss it? 2) Is there any implementation you know off that is suitable for wrapping or a good reference for an algorithm (I can find some, but if there's something well known...) Yours, Yosef.

4 6

module docstrings
by Ralf Gommers Nov. 4, 2009

Nov. 4, 2009

Hi all, The doc standard does not say anything about module docstrings, and before starting to clean those up it would be nice to agree on a format. Right now there are modules that: - do not have any docstring (matrixlib, fromnumeric, ..) - have license/author info in the docstring (scipy.interpolate.fitpack, ..) - list their routines (linalg, core, lib, ..) * some of these have a routines rst doc: routines.linalg.rst exists, routines.lib.rst does not. - have a single summary line (distutils, ..) - have some basic explanation (ma, scimath, ..) What do you all think about the following? : 1. every module needs a docstring with at least a summary line. 2. follow the doc standard where relevant, so typically: - summary line - extended summary paragraph(s) - see also, if needed - notes, if needed - references, if needed 3. no routine listing in the docstring, each module gets a corresponding routines.module.rst file 4. license and author info can be in the source file, but not in the docstring. (same principle for author info in reST docs by the way). Cheers, Ralf

5 9

scipy.test() failures
by Neil Killeen Oct. 29, 2009

Oct. 29, 2009

Hi all Hi - I am operating on an SGI Altix 3700 BX2 Itanium IA64 system running SUSE/Novelle Linux SLES10 - I am trying to install scipy V 0.7.1 - I am using python 2.4 - blas and lapack are installed via system rpms in /usr/lib soma:/usr/lib # ls -al *blas* *lapack* -rw-r--r-- 1 root root 929506 2007-04-14 02:16 libblas.a lrwxrwxrwx 1 root root 12 2008-03-12 12:20 libblas.so -> libblas.so.3 lrwxrwxrwx 1 root root 16 2008-03-12 12:20 libblas.so.3 -> libblas.so.3.0.0 -rwxr-xr-x 1 root root 678024 2007-04-14 02:16 libblas.so.3.0.0 lrwxrwxrwx 1 root root 20 2008-12-02 22:52 libgslcblas.so.0 -> libgslcblas.so.0.0.0 -rwxr-xr-x 1 root root 908816 2007-04-14 02:30 libgslcblas.so.0.0.0 -rw-r--r-- 1 root root 11208466 2007-04-14 02:16 liblapack.a lrwxrwxrwx 1 root root 14 2008-03-12 12:20 liblapack.so -> liblapack.so.3 lrwxrwxrwx 1 root root 18 2008-03-12 12:20 liblapack.so.3 -> liblapack.so.3.0.0 -rwxr-xr-x 1 root root 7503568 2007-04-14 02:16 liblapack.so.3.0.0 1. I did a vanilla build (python setup.py build/install) which appeared to work ok. 2. When I tried to run scipy.test() I saw errors. They are in two categories: ERROR: Failure: ImportError (/usr/local/lib/python2.4/site-packages/ scipy/linalg/fblas.so: undefined symbol: srotmg_) and ImportError: cannot import name {linalg, calc_lwork} Perhaps the latter are due to the former ? 3. I read in the FAQ that the solution to the srotmg error is to install blas myself either from netlib or via ATLAS 4. I decided to do the easy thing and installed blas from netlib in / usr/local/lib/libblas.a 5. I then deleted the scipy installation (/usr/local/lib/python2.4/ site-packges/scipy) and rebuilt and installed afresh. In the build log one can see that it appears to have found the /usr/ local/lib/libblas.a library E.g. blas_info: FOUND: libraries = ['blas'] library_dirs = ['/usr/local/lib'] language = f77 FOUND: libraries = ['blas'] library_dirs = ['/usr/local/lib'] define_macros = [('NO_ATLAS_INFO', 1)] language = f77 6. Then I - set env variable BLAS=/usr/local/lib/libblas.a (not sure useful since static library so should be no run-time dependency) - ran the test again with the same errors resulting as the initial. Can someone advise what more I must do to eradicate these run-time errors ? thanks Neil p.s. tried to post this last week but message did not arrive (perhaps as I included attachments which I have left off this time)

4 9

L0 norm
by Jason Rennie Oct. 26, 2009

Oct. 26, 2009

SciPy.linalg.basic.norm (version 0.7.0) raises a division-by-zero error when an ord value of zero is used. Yet, the L0 norm is commonly defined as the number of non-zero elements in a vector (see wikipedia link below). Has there been any discussion of adding an if statement to scipy.linalg.basic.norm to handle this case? Or, is this the behavior in a future version of scipy? If this functionality hasn't yet been added, would anyone object to a ticket being created requesting this feature? http://en.wikipedia.org/wiki/L0_norm#Zero_norm Cheers, Jason -- Jason Rennie Research Scientist, ITA Software 617-714-2645 http://www.itasoftware.com/

2 2

RFR: Proposed fixes in scipy.stats functions for calculation of variance/error/etc.
by Ariel Rokem Oct. 26, 2009

Oct. 26, 2009

Hi everyone, I have been working on some fixes to the functions in scipy.stats which calculate variance/error and related quantities. In particular, in order to comply with the deprecation warnings that appear in use of scipy.stats.samplevar/scipy.stats.samplestd, I have replaced use of these functions with calls to np.std/np.var. I have also cleaned up the documentation a bit. This can all be found here: http://codereview.appspot.com/141051 Cheers, Ariel -- Ariel Rokem Helen Wills Neuroscience Institute University of California, Berkeley http://argentum.ucbso.berkeley.edu/ariel

4 14

chararray docstrings
by Michael Droettboom Oct. 26, 2009

Oct. 26, 2009

I was able to make my big chararray commit today. If I understand correctly, I need to wait 24 hours for the doc editor to sync with SVN, and then I should mark all the chararray-related docstrings as "needs review". The primary change to the docstrings is that all of the methods of the chararray class are now free functions. These free functions represent the "primary" entry points, and thus have detailed documentation, and the chararray methods now have short "pointer" docstrings to the free functions. Where the docstring content itself has been updated, it is mainly to bring them closer to the Python standard library descriptions of these functions, which in most cases was more precise (since we are, in fact, calling the stdlib function under the hood) and concise (because the stdlib docs have been through a number of revisions and really get it right by now). I do have a concern about one phrase that was used in a number of places that probably deserves some discussion: "The chararray module exists for backwards compatibility with Numarray, it is not recommended for new development. If one needs arrays of strings, use arrays of dtype <http://docs.scipy.org/numpy/docs/numpy.dtype/#dtype> object." There are many use cases (such as handling a binary structured format like FITS) where a dtype of 'string_' is more appropriate than a dtype of 'object_', and we shouldn't imply that all uses of chararray should now use object arrays. Additionally, fast vectorized string operations will perform best on arrays of type 'string_' and 'unicode_', though 'object_' will work, it requires casting all objects to strings along the way, and could fail thousands of items in to an operation. It's a "best tool for the job" judgment call, not a "one tool fits all". Perhaps the above should read: "If one needs arrays of strings, use arrays of dtype <http://docs.scipy.org/numpy/docs/numpy.dtype/#dtype> string_ or unicode_. If one needs arrays of variable-length strings, use arrays of dtype object_." Mike -- Michael Droettboom Science Software Branch Operations and Engineering Division Space Telescope Science Institute Operated by AURA for NASA

6 21