Mailman 3 November 2010 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck 25 Nov '13

25 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

Subversion scipy.stats irregular problem with source code example
by James Phillips 10 Dec '10

10 Dec '10

Since I observed the following behavior in the SVN respository version of SciPy, it seemed to me proper to post to the dev mailing list. I'm using Ubuntu Lucid Lynx 32 bit and a fresh GIT of Numpy. I am not sure if a Trac bug report needs to be entered. Below is some example code for fitting two statistical distributions. Sometimes the numpy-generated data is fit, as I can see the estimated and fitted parameters are different. Sometimes I receive many messages repeated on the command line like: Warning: invalid value encountered in absolute Warning: invalid value encountered in subtract and the estimated parameters equal the fitted parameter values, indicating no fitting took place. Sometimes I receive on the command line: Traceback (most recent call last): File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1987, in func sk = 2*(b-a)*math.sqrt(a + b + 1) / (a + b + 2) / math.sqrt(a*b) ValueError: math domain error Traceback (most recent call last): File "example.py", line 10, in <module> fitStart_beta = scipy.stats.beta._fitstart(data) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/stats/distributions.py", line 1992, in _fitstart a, b = optimize.fsolve(func, (1.0, 1.0)) File "/home/zunzun/local/lib/python2.6/site-packages/scipy/optimize/minpack.py", line 125, in fsolve maxfev, ml, mu, epsfcn, factor, diag) minpack.error: Error occured while calling the Python function named func and program flow is stopped. In summary, three behaviors: (1) Fits OK (2) Many exceptions with no fitting (3) minpack error. Running the program 10 times or so will reproduce these behaviors without fail from the "bleeding-edge" repository code. James Phillips ######################################################## import numpy, scipy, scipy.stats # test uniform distribution fitting data = numpy.random.uniform(2.0, 3.0, size=100) fitStart_uniform = scipy.stats.uniform._fitstart(data) fittedParameters_uniform = scipy.stats.uniform.fit(data) # test beta distribution fitting data = numpy.random.beta(2.0, 3.0, size=100) fitStart_beta = scipy.stats.beta._fitstart(data) fittedParameters_beta = scipy.stats.beta.fit(data) print print 'uniform._fitstart returns', fitStart_uniform print 'fitted parameters for uniform =', fittedParameters_uniform print print 'beta._fitstart returns', fitStart_beta print 'fitted parameters for beta =', fittedParameters_beta print

7 42

ARPACK fixes before 0.9
by Pauli Virtanen 05 Dec '10

05 Dec '10

Hi, The Arpack interface in Scipy needs some fixes, which would be useful to get in for Scipy 0.9. I have them in a branch here: https://github.com/pv/scipy-work/tree/bug/1313-arpack Comments would be appreciated (esp. from David who wrote the original interface). I'm going to merge them probably the next W/E unless objections arise. Bugfixes: - Raise ArpackNoConvergence when ARPACK iteration does not converge. Previously, there was no way to know if convergence was obtained. - Fix a bug in return value extraction from dneupd, which resulted to invalid eigenvectors/eigenvalues being returned on non-convergence. - Allow complex matrices in the sparse SVD. It's a naive approach, but probably still better than nothing. Other changes: - Rename routines eigen* -> eigs* to avoid shadowing the "eigen" module. "eigs" is also the name for the equivalent routines in Octave/Matlab. - Remove the ``speigs`` ARPACK interface. It does not seem to make sense two have two different interfaces to the same library. - Un-deprecate .dot() method of sparse matrixes -- ndarrays also have it, so it makes sense to retain it. -- Pauli Virtanen

3 5

Sporadic failures of tests of signal.correlate with dtype complex64
by Warren Weckesser 30 Nov '10

30 Nov '10

On Mac OSX 10.5.8, I'm seeing occasional failures like the following: $ python -c "import scipy.signal; scipy.signal.test()" Running unit tests for scipy.signal NumPy version 1.5.0.dev8716 NumPy is installed in /Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy SciPy version 0.9.0.dev6856 SciPy is installed in /Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy Python version 2.6.5 |EPD 6.2-2 (32-bit)| (r265:79063, May 28 2010, 15:13:03) [GCC 4.0.1 (Apple Inc. build 5488)] nose version 0.11.3 ................./Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy/signal/filter_design.py:256: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless "results may be meaningless", BadCoefficients) ......................................................F............................................................................................................................................................................................................................................... ====================================================================== FAIL: test_rank1_same (test_signaltools.TestCorrelateComplex64) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy/signal/tests/test_signaltools.py", line 606, in test_rank1_same assert_array_almost_equal(y, y_r) File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy/testing/utils.py", line 774, in assert_array_almost_equal header='Arrays are not almost equal') File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy/testing/utils.py", line 618, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not almost equal (mismatch 10.0%) x: array([-6.76370811-8.55324841j, 0.68672836-4.2681613j , -3.22760987-8.69287109j, 0.75051951-5.50820398j, -7.33016682-1.14685655j, -5.99573374+7.84123898j,... y: array([-6.76370859-8.55324745j, 0.68672895-4.2681613j , -3.22761011-8.69286919j, 0.75051963-5.50820446j, -7.33016682-1.14685678j, -5.99573517+7.84123898j,... ---------------------------------------------------------------------- Ran 311 tests in 2.307s FAILED (failures=1) $ python -c "import scipy.signal; scipy.signal.test()" Running unit tests for scipy.signal NumPy version 1.5.0.dev8716 NumPy is installed in /Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy SciPy version 0.9.0.dev6856 SciPy is installed in /Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy Python version 2.6.5 |EPD 6.2-2 (32-bit)| (r265:79063, May 28 2010, 15:13:03) [GCC 4.0.1 (Apple Inc. build 5488)] nose version 0.11.3 ................./Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy/signal/filter_design.py:256: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless "results may be meaningless", BadCoefficients) .......................................................F.............................................................................................................................................................................................................................................. ====================================================================== FAIL: test_rank1_same_old (test_signaltools.TestCorrelateComplex64) ---------------------------------------------------------------------- Traceback (most recent call last): File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy/testing/decorators.py", line 257, in _deprecated_imp f(*args, **kwargs) File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/scipy/signal/tests/test_signaltools.py", line 641, in test_rank1_same_old assert_array_almost_equal(y, y_r) File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy/testing/utils.py", line 774, in assert_array_almost_equal header='Arrays are not almost equal') File "/Users/warren/tmp_py_install/lib/python2.6/site-packages/numpy/testing/utils.py", line 618, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not almost equal (mismatch 10.0%) x: array([ 2.46665049+2.02072477j, -7.42591763-0.54789257j, 3.41454220-0.15863085j, -0.14030695+5.01129198j, -2.11230707+2.68583822j, 7.78784609+7.19434834j,... y: array([ 2.46665049+2.02072501j, -7.42591763-0.54789257j, 3.41454196-0.15863061j, -0.14030659+5.01129246j, -2.11230707+2.68583822j, 7.78784752+7.19434786j,... ---------------------------------------------------------------------- Ran 311 tests in 2.623s FAILED (failures=1) The above tests are part of a suite of tests that use random data, and usually the tests all pass. It took several tries to get the above failures. I suspect the problem is simply that the default tolerance of 'assert_array_almost_equal' is too small for the complex64 data type for these tests. Could someone verify that they can reproduce those failures? Does simply increasing the tolerance of the test look like a reasonable fix? Warren

5 5

class.rst information and bug with the newer sphinx?
by Skipper Seabold 29 Nov '10

29 Nov '10

I've been playing around a bit with numpy/scipy's doc/source/_templates/autosummary/class.rst. First does this still work as expected with the newer sphinx? See for instance the method of the dtype class for numpy 1.4 (assumed built with older sphinx before 1.0?) and numpy 1.5 http://docs.scipy.org/doc/numpy-1.4.x/reference/generated/numpy.dtype.html#… http://docs.scipy.org/doc/numpy-1.5.x/reference/generated/numpy.dtype.html#… The newbyteorder method does not have a summary in the 1.5 version. Second, Is there anywhere to get more information on the autosummary templates? This covers layout.html <http://sphinx.pocoo.org/templating.html#templating> This mentions the possible files <http://sphinx.pocoo.org/ext/autosummary.html?highlight=templating#customizi…> but I am missing the details. In <http://svn.scipy.org/svn/numpy/trunk/doc/source/_templates/autosummary/clas…> where do the lines {% for item in methods %} {{ name }}.{{ item }} come from? Specifically, where does 'name' come from? Is this from Python, Sphinx, Numpy extension, or Jinja? Any help to demystify this would be appreciated. Skipper

2 2

deprecations/removals for 0.9.0
by Ralf Gommers 28 Nov '10

28 Nov '10

Hi all, Here is a reminder that if you want to deprecate things for 0.9.0, to do that in the next two weeks. Furthermore, I took a look at currently deprecated things that should be removed. Below is a long list, I believe that all items on it have been deprecated long enough that they can be removed. Please check and if something can *not* be removed, say so in the next few days. Also, if you deprecate more things, please state in the source code (as well as in the release notes) at what version removal should take place. Which reminds me: what should the version number after 0.9.0 be? Cheers, Ralf REMOVALS Stats ----- pdf_moments pdfapprox z zs stderr samplestd std var samplevar mean median cov corrcoef erfc Sparse ------ spkron speye spidentity dims and nzmax keywords in _sc_matrix constructor Sparse matrix methods: getdata rowcol listprint matvec matmat dot rmatvec Ndimage ------- output_type keyword in interpolate.py Linalg ------ the econ keyword in qr Signal ------ a whole bunch of stuff that I assume Warren is going to clean.

3 7

mstats.linregress return values (#1273)
by Ralf Gommers 28 Nov '10

28 Nov '10

Hi, stats.linregress returns five values: slope, intercept, r_value, p_value, sterr. mstats.linregress should return the same according to its docstring (and my expectation), but adds a sixth return value: Syy/Sxx, not clear what for. We can either document the return value, or just remove it. I'd think the latter makes more sense. Ralf

2 2

Fwd: [SciPy-User] FIR filter with arbitrary frequency response
by Warren Weckesser 25 Nov '10

25 Nov '10

Forwarded to the mailing list (didn't "reply to all" the first time). ---------- Forwarded message ---------- From: Warren Weckesser <warren.weckesser(a)enthought.com> Date: Sun, Nov 21, 2010 at 5:00 PM Subject: Re: [SciPy-User] FIR filter with arbitrary frequency response To: Nathaniel Smith <njs(a)pobox.com> Nathanial, thanks for the great feedback. Comments below... On Sat, Nov 20, 2010 at 11:25 PM, Nathaniel Smith <njs(a)pobox.com> wrote: > [Moving this to scipy-dev] > > On Sat, Nov 20, 2010 at 8:09 PM, Warren Weckesser > <warren.weckesser(a)enthought.com> wrote: > > There is one implemented as the function firwin2 currently under review > in > > this ticket: > > http://projects.scipy.org/scipy/ticket/457 > > It will eventually be added to scipy.signal. Included in the ticket is a > > patch file that can be applied to the latest version of the scipy source, > > and also a copy of just the updated file fir_filter_design.py. You could > > grab that and try it stand-alone (but you would have to comment out the > > local import of sigtools, which is used by the remez function in > > fir_filter_design.py). > > > > Feedback would be appreciated, so if you try it, be sure to write back > with > > comments or questions. > > Ha, what luck! Thanks! > > Some quick thoughts just based on reading the source code: > > -- Your handling of discontinuities in the desired frequency response > looks questionable to me (not that I'm an expert). Octave's version of > this function has a 'ramp_n' parameter that controls how much they > smooth out discontinuities, documented as "transition width for jumps > in filter response. Defaults to grid_n/20; a wider ramp gives wider > transitions but has better stopband characteristics." Your function, > OTOH, always uses the narrowest possible ramp. > > Octave's choice to make ramp_n defined in number-of-grid-point units > seems bizarre to me; I wouldn't copy that detail. But it probably > would be good to have some way to specify smarter handling of > discontinuities? > > As you noticed, currently it does not attempt to smooth the requested frequency response. If the user wants a ramp, they can just define their function with a ramp instead of a discontinuity. :) To handle transition width issues, I've used the Kaiser design method, which provides formulas for the number of taps and the Kaiser parameter beta in terms of the desired attenuation and transition width. This works reasonably well. I have nothing against adding some sort of "auto-smoothing" option. I'd prefer something with a sound mathematical basis, rather than an ad hoc tweak. But if there is a widely used heuristic that works, it should probably be included in the function. In the Octave example, what happens if the user includes a ramp in the response that isn't as wide as the transition width defined by grid_n/20? > -- I suspect the first part of the docstring: > > " From the given set of frequencies and gains, the desired response is > constructed in the frequency domain. The inverse FFT is applied to the > desired response to create the associated convolution kernel, and the > first `numtaps` coefficients of this kernel, scaled by `window`, are > returned." > > will not make sense to anyone who doesn't already know how this works. > > If it were me, I'd start by saying: > "Given any desired combination of frequencies (`freq`) and gains > (`gain`), this function constructs an FIR filter with linear phase and > (approximately) the given frequency response." > And I'd put the terse description of how it accomplishes this trick > lower down in the docstring. > > Thanks--after rereading, I see that it does need to be improved. > - In info.py, I'd expand the description of firwin to "Windowed FIR > filter design, with standard high/low/band-pass frequency response", > to make the contrast between it and firwin2 clear. > > Good idea. > - I'd also add a reference to http://www.dspguide.com/ch17/1.htm; > obviously any good DSP book will explain this idea, and the reference > you already have is fine, but it's much easier to click a link than to > go to the library... and the explanation is much less terse than the > one in the docstring :-). > > I'm OK with web links, but the numpy/scipy documentation guidelines state "Referencing sources of a temporary nature, like web pages, is discouraged." > -- Your doc string still refers to this function as 'fir2' at one point > > Thanks, I'll fix that. > -- I have a mild preference for specifying the sampling frequency > instead of the nyquist frequency, but that's not a big deal. What's > really bad is that we have no standard for this (it looks like the > only other filter design function that lets you work in sampling units > instead of normalized units is remez, and it uses the terribly named > 'Hz'). Maybe I should start another thread to discuss that... > > Personally, I'm fine with either sample rate or Nyquist rate, but there are other functions where the default is to express frequencies relative to the Nyquist rate (e.g. all the IIR functions). I also noticed that remez has the unfortunate keyword 'Hz'. For the sake of consistency, I would like to deprecate that keyword, and add 'nyq' (or perhaps 'nyquist', or...) as the preferred keyword. This option should also be added to firwin(), so all these FIR filter design functions have a consistent interface. > I don't have any authority or anything, but with those issues fixed > I'd vote for committing it. Tests seem reasonably complete. > > Great--thanks again for the comments. Warren

2 2

incorrect function-call count in fsolve
by Andrew Hawryluk 24 Nov '10

24 Nov '10

I have just noticed that scipy.optimize.fsolve reports two fewer function calls than actually occur. It returns the number of Fortran function calls, but the Python code calls the function once to check the dimensions and the C code calls the function again for error checking (assuming I have understood the code correctly). The same thing is true for calls to the Jacobian. My two questions are: Would it be appropriate for fsolve to return the number of Fortran calls + 2? Would it be appropriate to provide an 'unsafe' option that bypasses these calls? For some of my work I am calling fsolve several times within other optimization/solver loops, so the Fortran code may only be needing 1 Jacobian and 2 or 3 function calls much of the time. Bypassing the additonal 2 Jacobians and function calls could speed things along a bit. Cheers, Andrew from scipy.optimize import fsolve import numpy as np nfev = 0 def func(x): global nfev nfev += 1 return np.cos(x) solution = fsolve(func, 1.3, full_output=1) print 'fsolve performed %d function calls' % solution[1]['nfev'] print 'but there were %d function calls in total' % nfev

1 0

Pythonic Simulink clone... - how much interest?
by burley 22 Nov '10

22 Nov '10

Visual modeling support as offered by the Simulink/Matlab and Xcos/Scilab suites is very powerful, for the modeler as well as in communicating concepts and ideas to others. Neither is Pythonic or extendable by Python, unfortunately (disregarding the prohibitively expensiveness of Sim/Matl). Wouldn't it be *very* nice to have a free (LGPL++?), user extensible, Python based Simulink alternative, based on Numpy, SciPy, Wx and/or Qt-like, and the rest. All (e.g. most of) the algorithms are there, many would benefit. I know of Viper (Tcl/Tk), Elefant (Bayesian Learning), ... but none of them is general purpose or flexible in its use. Only a collaborative effort may make that happen I would believe. What's the community interest in such an effort? (I for myself am not a GUI programmer...). -- Burley

8 12