Mailman 3 April 2008 - 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

Question about Kaiser Implementation in firwin
by Buehler, Eric (AGRE) 04 Oct '08

04 Oct '08

Hello, I have been using the signal.firwin FIR design tool and I have a question about the implementation. Scipy-0.4.8 Numpy-0.9.6 I am generating a kaiser window with the following parameters: N = 128 Cutoff = ~4e7 width = .3 All of the special kaiser generation parameters in the function follow Oppenheim and Schafer well. Even the sinc function for linear phase is follows exactly. However, the last line of the firwin function is causing me some heartburn. filter_design.py 1538 win = get_window(window,N,fftbins=1) 1539 alpha = N//2 1540 m = numpy.arange(0,N) 1541 h = win*special.sinc(cutoff*(m-alpha)) 1542 return h / sum(h) Line 1542 of filter_design.py, "return h / sum(h)", normalizes the function where it doesn't seem necessary, at least in the kaiser window case. Without the normalization, the kaiser window already returns a value of 1 at the zero frequency point. This normalization scales all of the data, making the window difficult to use in the frequency domain. Can someone point me to the rationale for this line? Looking at the code, this seems to be a pretty recent change (within the last year/year and a half). Thanks, Eric Buehler eric <dot> buehler <at> smiths-aerospace <dot> com ****************************************** The information contained in, or attached to, this e-mail, may contain confidential information and is intended solely for the use of the individual or entity to whom they are addressed and may be subject to legal privilege. If you have received this e-mail in error you should notify the sender immediately by reply e-mail, delete the message from your system and notify your system manager. Please do not copy it for any purpose, or disclose its contents to any other person. The views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of the company. The recipient should check this e-mail and any attachments for the presence of viruses. The company accepts no liability for any damage caused, directly or indirectly, by any virus transmitted in this email. ******************************************

3 2

[ANN] EuroSciPy Abstracts Deadline Reminder
by Travis Vaught 30 Apr '08

30 Apr '08

Greetings, Just a reminder: the abstracts for the EuroSciPy Conference in Leipzig are due by midnight tonight (CST, US [UTC -6]) April, 30. If you'd like to present, please submit your abstract as a PDF, MS Word or plain text file to euroabstracts(a)scipy.org. For more information on the EuroSciPy Conference, please see: http://www.scipy.org/EuroSciPy2008

1 0

random seed in scipy.stats
by Hoyt Koepke 29 Apr '08

29 Apr '08

Hello, I may have missed something here, so I'll first ask: Is there a way to specify a random number seed in scipy.stats? I'm building a project using it which requires exactly reproducible distribution samples, so being able to specify a random seed is a must. I've been looking around in the code, and I don't see a way to do it. Even if it is there, it calls numpy.random directly, which would not work reliably. However, it seems that it'd be easy to add in by initializing a RandomState object in the constructors of both rv_continuous and rv_discrete, seeding it using an optional kw arg that defaults to None. All the code that currently draws samples from numpy.random directly would then draw samples from this. As far as I can see, this wouldn't break backwards compatibility and would add in a useful feature. The only change (I think) would be to distributions.py. Have I missed anything? If this is acceptable, I'll gladly make the change and post a diff and a unit test. --Hoyt -- +++++++++++++++++++++++++++++++++++ Hoyt Koepke UBC Department of Computer Science http://www.cs.ubc.ca/~hoytak/ hoytak(a)gmail.com +++++++++++++++++++++++++++++++++++

2 10

Changeset 4181
by Nils Wagner 28 Apr '08

28 Apr '08

Hi All, Changeset 4181 introduces ====================================================================== ERROR: test_construction (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 168, in test_construction P = PiecewisePolynomial(self.xi,self.yi,3) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_derivative (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 175, in test_derivative P = PiecewisePolynomial(self.xi,self.yi,3) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_derivatives (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 178, in test_derivatives P = PiecewisePolynomial(self.xi,self.yi,3) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_incremental (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 199, in test_incremental P.append(self.xi[i],self.yi[i],3) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_scalar (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 171, in test_scalar P = PiecewisePolynomial(self.xi,self.yi,3) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_shapes_scalarvalue (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 203, in test_shapes_scalarvalue P = PiecewisePolynomial(self.xi,self.yi,4) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_shapes_scalarvalue_derivative (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 209, in test_shapes_scalarvalue_derivative P = PiecewisePolynomial(self.xi,self.yi,4) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar ====================================================================== ERROR: test_vector (test_polyint.CheckPiecewise) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/tests/test_polyint.py", line 189, in test_vector for i in xrange(len(ys[0][0]))] File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 448, in __init__ self.extend(xi[1:],yi[1:],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 539, in extend self.append(xi[i],yi[i],orders) File "/usr/local/lib64/python2.5/site-packages/scipy/interpolate/polyint.py", line 503, in append raise ValueError, "Each derivative must be a scalar" ValueError: Each derivative must be a scalar Nils

2 2

MIT license
by Nathan Bell 28 Apr '08

28 Apr '08

We've received some MIT-licensed code to read matrices in the Harwell-Boeing file format[1]. My understanding is that the MIT license is essentially the same as the BSD, so is there any need to relicense the code? Generally, can MIT licensed code be included in SciPy without further consideration? [1] http://projects.scipy.org/scipy/scipy/ticket/354 Dominique, thanks for your contribution and sorry for the slow progress on this ticket. -- Nathan Bell wnbell(a)gmail.com http://graphics.cs.uiuc.edu/~wnbell/

5 4

Re: [SciPy-dev] Status of sandbox.spline and my other code.
by John Travers 26 Apr '08

26 Apr '08

---------- Forwarded message ---------- From: John Travers <jtravs(a)gmail.com> Date: Tue, Apr 8, 2008 at 9:55 AM Subject: Re: [SciPy-dev] Status of sandbox.spline and my other code. To: Jarrod Millman <millman(a)berkeley.edu> On Mon, Apr 7, 2008 at 11:56 AM, Jarrod Millman <millman(a)berkeley.edu> wrote: > > Hey John, > > I am hosting a NumPy/SciPy sprint this week at UC Berkeley and among > other things am hoping to make some progress on cleaning up the > various interpolation/spline code in SciPy as well as continue > removing the scipy sandbox. Do you have any more plans for your > sandbox.spline code? Any thoughts on where it should go or what still > needs to be done? > The code in sandbox.spline should be a complete replacement for the code in scipy.interpolate.fitpack and fitpack2 including the f2py wrappers of the fitpack fortran code (i.e. replacing the hand crafted c wrappers). Apart from using only f2py, a number of improvements were made to the docstrings and error checking and some bugs fixed and quite a few unit tests added. Originally I was advocating that we have a separate spline package with the fitpack stuff and a separate interpolate package that can utilise either the spline or rbf of delaunay etc. packages as these can all be used for more than just interpolation. This didn't seem to gain any traction on the mailing list, so I left it. I would still advocate this, but recognise the reasons that people might not want to expand the number of packages so much. I left this code about a year ago and then got bogged down with my thesis work which I finished a few weeks back. If they code/improvements might be used then I can find the time this week to check it is all working and still up to date. Then I could either supply a set of patches against the main interpolate code which could be more easily understood, or provide a mercurial/bzr branch for someone else to merge from (if the bzr/mercurial stuff is setup yet). Let me know if that might be useful. There are also further things that were on my todo list for interpolate, including using RectBivariateSPline where possible as it is much better for interpolation than the currently used surfit. (I have a unit test which segfaults the current code, but not teh RectBivariateSpline code). Hope this helps, John

3 3

2D array to 3D
by tournesol 26 Apr '08

26 Apr '08

Hi All. Is there a easy way to insert 1D(j) array into another 2D array(B:jxk) and conver B to B:ixjxk ? ex:) >>> from numpy import * >>> a=arange(4) >>> a array([0, 1, 2, 3]) >>> b=arange(9) >>> b.shape=3,3 >>> b array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) I just wanna insert A into B B:1x3x3, [[[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]] B:2x3x3, [[[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]] 　　 [[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]]] B:3x3x3, [[[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]] [[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]] [[ 0, 1, 2, 3], [ 4, 5, 6, 7], [ 8, 9, 10, 11]]] Thanks for your help

1 0

Interesting ARPACK branch
by Neilen Marais 25 Apr '08

25 Apr '08

Hi, I've found that figuring out what is the latest ARPACK release and weather you've got all the newest patches installed is quite a mission. It seems someone has take up a little of the the maintainership issue: http://mathema.tician.de/software/arpack Might be useful to pull the ARPACK code for future scipy releases from there. Just my lazy web-procrastination-induced 2c Neilen

2 1

Inclusion of cython code in scipy
by Anne Archibald 25 Apr '08

25 Apr '08

Hi, What do people think about including cython code in scipy? Pros: * Easy way to write faster code * Convenient way to interface to compiled code Cons: * Yet another interfacing setup * Yet another language readers of the code need to understand * Additional dependency for developers * cython is in flux and drastic improvements in its handling of numpy arrays are hoped for * Requires support in build tools (distutils/makefiles/scons magic) For one data point, a quick and nasty cythonization of my polynomial interpolation code got something like a 4x speedup for large arrays of evaluation points. This is a big improvement but far from the 20-100x cython claims for raw python. I can't judge the behaviour for shorter arrays without doing a more reasonable job of converting to cython, but by looking at the (horrible) generated C code, the inner loops are raw C with no python objects at all. (The 4x figure also uses gcc's -march and -O options.) I take this as evidence that numpy's arrays are doing a fairly good job of handing off the heavy lifting to compiled code. I can definitely see that if a little work went into cython's support of numpy arrays specifically the code could get a lot cleaner: things like c[i*k+j] could become c[i,k] with support for arbitrary strides and (possibly) type conversion. So maybe it doesn't make sense to start producing lots of cython code until that happens. Anne

19 54