Mailman 3 September 2011 - 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

prepending underscores to private module names
by Ralf Gommers 16 May '12

16 May '12

Hi, A while ago we discussed making a clear distinction between public and private modules by using underscores consistently, see http://thread.gmane.org/gmane.comp.python.scientific.devel/14936. I've done this for a few modules now, see https://github.com/rgommers/scipy/tree/underscores. Most files could simply be underscored, some other would conflict with extension modules so I also had to append _py to the name. Also I added a file doc/source/api.rst that states which modules are public (up for discussion of course). With google code search I checked how common it is to import from modules that are supposed to be private, for example "from scipy.optimize.optimize import fmin". It's not very common but does happen, so deprecation warnings should be put in. Before going further I'd like to check that this looks okay to people. What do you think? Cheers, Ralf <https://github.com/rgommers/scipy/tree/underscores>

4 4

studentized range approximations
by Roger Lew 28 Mar '12

28 Mar '12

Hi, I have implemented some approximations for studentized range quantiles and probabilities based on John R. Gleason's (1999) "An accurate, non-iterative approximation for studentized range quantiles." Computational Statistics & Data Analysis, (31), 147-158. Probability approximations rely on scipy.optimize.fminbound. The functions accept both scalars or array-like data thanks to numpy.vectorize. A fair amount of validation and testing has been conducted on the code. More details can be found here: http://code.google.com/p/qsturng-py/ I welcome any thoughts as to whether you all think this might be useful to add to SciPy or make into a scikit. Any general comments would be helpful as well. I should mention I'm a cognitive neuroscientist by trade, my use of statistical jargon probably isn't that good. Regards, Roger Roger Lew

3 5

installation of scipy on Mac OS X 10.7
by Ben Willmore 26 Jan '12

26 Jan '12

I have spent a day attempting to get a version of scipy that passes unit tests on Mac OS X 10.7. I hope the results of my investigations (below) are useful to others trying to do the same thing. Ben == Introduction I have attempted to obtain a version of scipy on Mac OS X 10.7 that passes unit tests, both by downloading it, or by compiling against the system's python, using the three C compilers available, and various compile flags. == Conclusions 1. Both Enthought Python 7.1 and the scipy superpack crash during scipy unit tests. 2. Compilation with the default C compiler (llvm-gcc-4.2) also leads to a crash during unit tests. 3. Compilation using either of the other C compilers (gcc-4.2 or clang) allows scipy to complete unit tests, only if the fortran flag --ff2c is used. In both cases, scipy passes almost all tests -- producing errors on one, and failing two more. With gcc-4.2, numpy passes all tests. With clang, numpy fails one test. 4. The following recipe was most successful for me: * Remove everything in /usr/local and /Library/Python/2.7/site-packages (BEWARE! These directories may well contain stuff you care about) * Install gfortran from here: http://r.research.att.com/gfortran-lion-5666-3.pkg mkdir ~/tmp cd ~/tmp git clone git://github.com/numpy/numpy.git cd numpy export CC=gcc-4.2 export CXX=g++-4.2 export FFLAGS=-ff2c python setupegg.py build --fcompiler=gfortran sudo python setupegg.py install cd ~/tmp git clone git://github.com/scipy/scipy.git cd scipy export CC=gcc-4.2 export CXX=g++-4.2 export FFLAGS=-ff2c python setupegg.py build --fcompiler=gfortran sudo python setupegg.py install (repeat for matplotlib, ipython ...) == Details of test results = Enthought python 7.1 numpy.test(): OK scipy.test(): Crash (Segmentation fault) = Scipy superpack numpy.test(): 3 failures: FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, 1, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, -1, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, 0.0, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, -0.0, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, inf, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, -inf, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, nan, inf, inf, inf) FAIL: test_umath_complex.TestCsqrt.test_special_values(<ufunc 'sqrt'>, -inf, 1, 0.0, inf) scipy.test(): Crash (Abort trap) = Compilation from scratch (see above for recipe) = llvm-gcc-4.2 = (no compiler flags) numpy.test(): OK scipy.test(): Crash (Segmentation fault: 11) = gcc-4.2 = CC=gcc-4.2 CXX=g++-4.2 numpy.test(): OK scipy.test(): Hang (during ARPACK tests) = clang = CC=clang CXX=clang++ numpy.test(): 1 failure: FAIL: Test basic arithmetic function errors AssertionError: Type <type 'numpy.complex64'> did not raise fpe error ''. scipy.test(): Hang (during ARPACK tests) = llvm-gcc-4.2 and -ff2c = FFLAGS=-ff2c numpy.test(): OK scipy.test(): Crash (Segmentation fault: 11) = gcc-4.2 and -ff2c = CC=gcc-4.2 CXX=g++-4.2 FFLAGS=-ff2c numpy.test(): OK scipy.test(): 1 error, 2 failures: ERROR: test_arpack.test_hermitian_modes(True, <gen-hermitian>, 'F', 2, 'SA', None, None, <function aslinearoperator at 0x1098cd500>) FAIL: test_iterative.test_convergence(<function bicg at 0x1098cdd70>, <nonsymposdef>) FAIL: test_expon (test_morestats.TestAnderson) = clang and -ff2c = CC=clang CXX=clang++ FFLAGS=-ff2c numpy.test(): 1 failure: FAIL: Test basic arithmetic function errors AssertionError: Type <type 'numpy.complex64'> did not raise fpe error ''. scipy.test(): 1 error, 2 failures: ERROR: test_arpack.test_hermitian_modes(True, <gen-hermitian>, 'F', 2, 'SA', None, None, <function aslinearoperator at 0x1098cd500>) FAIL: test_iterative.test_convergence(<function bicg at 0x1098cdd70>, <nonsymposdef>) FAIL: test_expon (test_morestats.TestAnderson)

6 6

On the leastsq/curve_fit method
by Gianfranco Durin 03 Oct '11

03 Oct '11

Dear all, I wanted to briefly show you the results of an analysis I did on the performances of the optimize.leastsq method for data fitting. I presented these results at the last Python in Physics workshop. You can download the pdf here: http://emma.inrim.it:8080/gdurin/talks. 1. The main concern is about the use of cov_x to estimate the error bar of the fitting parameters. In the docs, it is set that "this matrix must be multiplied by the residual standard deviation to get the covariance of the parameter estimates -- see curve_fits."" Unfortunately, this is not correct, or better it is only partially correct. It is correct if there are no error bars of the input data (the sigma of curve_fit is None). But if provided, they are used as "as weights in least-squares problem" (curve_fit doc), and cov_x gives directly the covariance of the parameter estimates (i.e. the diagonal terms are the errors in the parameters). See for instance here: http://www.gnu.org/s/gsl/manual/html_node/Computing-the-covariance-matrix-o…. This means that not only the doc needs fixing, but also the curve_fit code, those estimation of the parameters' error is INDEPENDENT of the values of the data errors in the case they are constant, which is clearly wrong. I have never provided a patch, but the fix should be quite simple, just please give me indication on how to do that. 2. The convergence of the fit in the most difficult cases (see page 15 of my presentation) can required up to about 3000 iterations, reduced to 800 when using analytical derivatives. For quite a long time I did not realized that the fit needed more iterations that the number set by maxfev, and thus I started to think that the leastsq was not good enough for 'hard' data. As a matter of fact, maxfev : int The maximum number of calls to the function. If zero, then 100*(N+1) is the maximum where N is the number of elements in x0. is pretty low, so I suggest to increase the prefactor 100 to 1000. A relatively huge number is not a problem, by the way, because if the system is sloppy. i.e. one parameter does not move too much, the routine stops and complains with "Both actual and predicted relative reductions in the sum of squares are at most 0.000000 and the relative error between two consecutive iterates is at most 0.000000" as in the case of boxBod at pg. 15. By the way, we should also advice that the in case of analytical derivative this number is half, even if I personally would keep the same number for both cases. This is all for the moment. Many thanks for your attention, and sorry for the long mail Gianfranco

5 11

Function for finding relative maxima of a 1d array
by Jacob Silterra 30 Sep '11

30 Sep '11

Hello all, I just opened a pull request to add scipy.signal._peak_finding.find_peaks, which finds the relative maxima in a 1d ndarray. The algorithm behind this function has been discussed on the scipy list, but here's a quick recap: The idea is basically convolve the function with a wavelet of varying widths (default is ricker, or "mexican hat" wavelet) and identify as peaks those locations where there are relative maxima (meaning they are larger than their nearest neighbor on each side) at each of several widths, and have a sufficiently high SNR. More details are provided in the docstrings of the functions. Pull request at https://github.com/scipy/scipy/pull/85. I believe the code quality/docstrings/test coverage are pretty good, however it would be of great benefit for someone with experience with similar algorithms to review the code. -Jacob PS I've sent this to both the scipy and numpy lists, apologies to those who get it twice.

1 0

lazy imports feedback request
by Paul Ivanov 30 Sep '11

30 Sep '11

Friends and uncles, I'm seeking feedback about some lazy import functionality I'm trying to implement in nitime. In short, "foo = LazyImport('foo')" returns something that will act as foo if you did just "import foo" - but only perform the import when foo is actually used, not when it is defined. I'll outline the highlights in this email, but you can see the full four-leaf clover pull request discussion here: https://github.com/nipy/nitime/pull/88 Here's the code, with an explanatory docstring class LazyImport(object): """ This class takes the module name as a parameter, and acts as a proxy for that module, importing it only when the module is used, but effectively acting as the module in every other way (including inside IPython with respect to introspection and tab completion) with the *exception* of reload(). >>> mlab = LazyImport('matplotlib.mlab') No import happens on the above line, until we do something like call an mlab method or try to do tab completion or introspection on mlab in IPython. >>> mlab <module 'matplotlib.mlab' will be lazily loaded> Now the LazyImport will do an actual import, and call the dist function of the imported module. >>> mlab.dist(1969,2011) 42.0 >>> mlab <module 'matplotlib.mlab' from '.../site-packages/matplotlib/mlab.pyc'> """ def __init__(self, modname): self.__lazyname__= modname def __getattribute__(self,x): # This method will be called only once name = object.__getattribute__(self,'__lazyname__') module =__import__(name, fromlist=name.split('.')) # Now that we've done the import, cutout the middleman class LoadedLazyImport(object): __getattribute__ = module.__getattribute__ __repr__ = module.__repr__ object.__setattr__(self,'__class__', LoadedLazyImport) return module.__getattribute__(x) def __repr__(self): return "<module '%s' will be lazily loaded>" %\ object.__getattribute__(self,'__lazyname__') The win is that when you get a faster import and smaller memory footprint by not loading those modules which your code path ends up not using. The functionality of lazyimports.LazyImport is generic enough to allow the lazily imported module to act as the module in almost every way (tab completion, introspection for docstrings and sources) except reloading is not supported. Although I have not spent a lot of time trying to make reloading work, I think it is not a big deal because here in the intent is to lazily load external packages (matplotlib, scipy, and numpy's nosetools) which are not meant to be developed simultaneously with the project wishing to load them in a lazy manner. A difference with another strategy, such as one Keith Goodman highlighted in a thread announcing Bottlneck 0.4.1 on [SciPy-User] back in March which looks like this: # http://wiki.python.org/moin/PythonSpeed/PerformanceTips#Import_Statement_Ov… email = None def parse_email(): global email if email is None: import email is that LazyImport keeps you from having to check if a given module has been lazily loaded in every place that it is used. Fernando Perez suggested I contact this list because he recalls similar functionality being implemented in scipy at some point, and then removed because it was problematic. thanks a lot in advance for any feedback. I'm cc-ing my primary address in hope that your mail reader will do the same. best, -- Paul Ivanov

2 1

Arpack problems and 0.10.0 release
by Ralf Gommers 25 Sep '11

25 Sep '11

Hi, The first release candidate for scipy 0.10.0 will have to be delayed a bit, due to problems with Arpack: http://projects.scipy.org/scipy/ticket/1515 (test failure on 64-bit Linux) http://projects.scipy.org/scipy/ticket/1523 (crash on 64-bit OS X) I'm hoping someone has the time and knowledge to fix this. For the OS X crash, it would also be helpful if someone could (a) reproduce it and (b) narrow it down to a self-contained example for easier debugging. Thanks, Ralf

1 0

Re: [SciPy-Dev] [SciPy-User] polylogarithm?
by Johann Cohen-Tanugi 25 Sep '11

25 Sep '11

For the record, I created a first pull request, with limited scope, in order to make the C call to cephes described below possible. Further details in https://github.com/scipy/scipy/pull/84 I will switch to the dev list from now on. best, Johann On 09/23/2011 06:18 PM, Pauli Virtanen wrote: > Fri, 23 Sep 2011 15:21:12 +0200, Johann Cohen-Tanugi wrote: > [clip] >> zetac(z) knows how to eat negative arguments, but not zeta(z,1).... Is >> there a reason why special.zeta does not default to 1+special.zetac for >> s=1? This would make the behavior of the 2 functions more identical. > Probably no reason, except that it wasn't implemented. mpmath is > impressive, and in several ways ahead of scipy.special --- or at least in > the parts where the problems overlap, as you can do tricks with arbitrary > precision that are not really feasible. > > Note that if you need to call the zeta function from the Cython extension, > call the C library directly: > > cdef extern from "cephes.h": > double zeta(double x, double q) > double zetac(double x) > > and link the extension with the "sc_cephes" library. > > *** > > There's a formula (look in the mpmath sources ;) for the transform > from x< 0 to x> 0 for zeta(x, a) for general a. But that needs polylog. > > The implementation for zetac(x) for x< 0 seems also a bit incomplete, > as it goes only down to -30.8148. It seems this is due to a silly reason, > it should use gammaln instead of Gamma to avoid the overflow. > > Pauli > > _______________________________________________ > SciPy-User mailing list > SciPy-User(a)scipy.org > http://mail.scipy.org/mailman/listinfo/scipy-user >

1 0

Docs.scipy.org wiki edit rights
by Gianfranco Durin 22 Sep '11

22 Sep '11

Hi all, I am extensively using the least-square fitting method and I recently found a bug in the documentation together with some improvements to add. I'd like to get the rights to edit the docs. I will soon present the results of my analysis to the list Many thanks Gianfranco

2 1