Mailman 3 December 2007 - 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

Moving numexpr to its own project
by David M. Cooke 15 Apr '08

15 Apr '08

Hi, With the pending closing of the scipy.sandbox (anybody want to cough up more details on that, btw? When? How?), I want to move numexpr to its own project, likely hosted on Google Code. It's not something I think is "sciencey" enough for a scikit (plus, I don't like sharing a bug database with other, separate, projects), and it's more like weave (but simpler). Tim H., Francesc, Ivan? You're the main contributors; does this work for you? -- |>|\/|< /------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm(a)physics.mcmaster.ca

6 12

Status of sandbox.spline and my other code.
by John Travers 07 Apr '08

07 Apr '08

Hi all, This email is intended to clarify the status of various bits of code I've committed/not committed so I don't leave too much of a mess after me. I'm about to start writing my PhD thesis, so I'll be too busy to work on scipy until September. 1. sandbox.spline stuff The module here was supposed to be a tidy up of scipy.interpolate. All of the dierkx functionality has been moved to f2py wrappers, which I think makes maintenance easier. However, it does not add any new functionality, and in retrospect appears to have been a waste of my time. So I'll leave to you to decide if you want to integrate it into scipy.interpolate or not; though it was originaly planned to be a seperate module to clear up the ambiguity between interpolation and smoothing spline functionality. However I have added 20 unit tests to the code (most of which simply check the wrapper against the pure fortran output, but still useful I think) which could easily be moved over to the current scipy.interpolate module. 2. Radial basis function module in sandbox (rbf) This code has recently had attention from Robert Hetland and is quite improved. I'm not sure if I will ever go into scipy though?? If not, I will put it into a scikit at some point. Related to this, the wiki page has been updated, however, it is at http://www.scipy.org/RadialBasisFunctions but should probably be in the Cookbook. I don't know how to move it over, so some pointers would be helpful. In addition there are four unused attachments I uploaded that could be deleted. 3. Boundary value ODE solver (BVP_SOLVER) I have fairly functional code interfacing to some very good BVP solver code (http://cs.smu.ca/~muir/BVP_SOLVER_Webpage.shtml). The license is good for scipy (as confirmed by the authors), but the code if fortran 95 - therefore cannot apparently go into scipy. I'll make a scikit when I get time (probably September). 4. An interface to the pikaia genetic algorithm routine I have a python interface to pikaia fortran code for genetic algorithm optimisation of real valued functions (http://www.hao.ucar.edu/Public/models/pikaia/pikaia.html). The fortran code itself is robust and highly optimised for numerical work, though it is nowhere near as general as the sandbox.ga module. I think it is useful to go in scipy.optimization alongside the annealing module. If people disagree I'll make it into a scikit. I hope this clears a few things up. Best regards, John Travers

4 4

ZVODE, a Double Precision Complex ODE Solver
by Nils Wagner 19 Feb '08

19 Feb '08

This might be of interest. Nils From: Alan Hindmarsh <alanh(a)llnl.gov> Date: Wed, 6 Dec 2006 18:20:32 -0500 Subject: ZVODE, a Double Precision Complex ODE Solver Of the many solvers freely available for ODE initial value problems, there seem to be very few, if any, for the case of double precision complex dependent variables, even though there have been occasional user calls for this functionality. (Of course a complex system can be treated as a real system of twice the size, but direct treatment as a complex system is preferable.) I have provided such a solver: ZVODE. ZVODE is a straightforward modification of the Fortran ODE solver DVODE, for double precision real systems, making it part of a trio, SVODE/DVODE/ZVODE. ZVODE is available in two locations: (1) Netlib, in /ode (see file zvode.f), and (2) the LLNL CASC software site, www.llnl.gov/CASC/software.html (see VODE). The latter site also includes a demonstration program, zvdemo. Alan Hindmarsh Center for Applied Scientific Computing Lawrence Livermore National Laboratory

4 5

The future of the scipy.sandbox and a reminder of upcoming doc-day
by Jarrod Millman 25 Jan '08

25 Jan '08

Hello, I started a blog, which has two posts about the NumPy/SciPy code sprint/strategic planning meeting early this month. It doesn't have very much content yet, but I thought it would be worth letting everyone know where it is: http://jarrodmillman.blogspot.com/ The most recent entry addresses some of the reasons we decided to retire the sandbox. I also created a ticket for removing the sandbox: http://projects.scipy.org/scipy/scipy/ticket/573 I will work on further documenting what went on at the sprint and some other things tomorrow during the first NumPy/SciPy doc-day. As a reminder, Friday's doc-day will be virtual where anyone with an internet connection can join in on the scipy channel on irc.freenode.net. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

7 25

Re: [SciPy-dev] Sparse BLAS in scipy
by Anand Patil 16 Jan '08

16 Jan '08

On Nov 28, 2007 4:06 PM, Nathan Bell <wnbell(a)gmail.com> wrote: > On Nov 27, 2007 4:35 PM, Anand Patil <anand.prabhakar.patil(a)gmail.com> > wrote: > > That and the analogous forward multiplication with a triangular matrix > are > > the two things I really need right now, yeah. The backend would be > simple, > > yet painful; upper/lower, transpose and side arguments, in addition to > > accounting for all the possible combinations of matrix formats, and the > > strides when 'B' is dense, would make for a huge number of bug > > opportunities. That's what made me think wrapping a library would be > less > > work than writing even 'two' routines from scratch. > > I guess I don't fully understand what you need here. What I imagined > was two methods for forward/backward solves with lower/upper > triangular matrices in CSR format. > I was thinking the triangular matrices could be any format, and the 'b' matrix/vector could be any format, in which case the multiplicity of functions rivals that of the ordinary multiplications involving sparse matrices. BTW, for such matrices, isn't this is equivalent to Gauss-Seidel[1] > sweep in the appropriate direction. If so, then we might just move > something like [2] to scipy.linsolve. > I can't see any difference... it would be nice to have the L3 version too, though, and in that case b could be in any format. Also, are you sure the sparse LU solvers wouldn't do this for you anyway? > Nope, will have to look more closely at those. > About a year ago I implemented sparsetools from scratch. I looked at > the same packages you listed before and came to the conclusion that > the amount of code needed was rather small and that by writing the > library to fit well with SciPy was worth the effort. > > Ultimately, I think making the backend code amenable to SWIG and > scipy.sparse saved more effort than starting with an existing library. OK, if you've done the research already and come to that conclusion I should just do it! Regardless of the backend (LU, gauss-seidel or from scratch) here's what I'm proposing for the Python interface: B=trimult(A, B, uplo='U', side='L', inplace=True) B=trisolve(A, B, uplo, side, inplace) A is triangular, B is rectangular or a vector. We could trim the uplo and side arguments if there were 'triangular matrix' classes that kept track of whether they're upper or lower, but that would be more trouble than it's worth, no? I'm thinking if inplace=True, B will only be overwritten if it's possible to recycle the memory efficiently, so the only safe usage is to use the return value. How does this sound? Anand

2 4

planet.scipy.org
by Jarrod Millman 31 Dec '07

31 Dec '07

Hey, I just wanted to announce that we now have a NumPy/SciPy blog aggregator thanks to Gaël Varoquaux: http://planet.scipy.org/ Feel free to contact me if you have a blog that you would like included. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

1 0

Re: [SciPy-dev] [Numpy-discussion] Doc-day
by Travis E. Oliphant 29 Dec '07

29 Dec '07

Stefan van der Walt wrote: > On Thu, Dec 27, 2007 at 09:27:09PM -0800, Jarrod Millman wrote: > >> On Dec 27, 2007 7:42 PM, Travis E. Oliphant <oliphant(a)enthought.com> wrote: >> >>> Doc-day will start tomorrow (in about 12 hours). It will be Friday for >>> much of America and be moving into Saturday for Europe and Asia. Join >>> in on the irc.freenode.net (channel scipy) to coordinate effort. I >>> imaging people will be in an out. I plan on being available in IRC from >>> about 9:30 am CST to 6:00 pm CST and then possibly later. >>> >>> If you are available at different times in different parts of the >>> world, jump in and pick something to work on. >>> >> Since this is our first doc-day, it will be fairly informal. Travis >> is going to be trying to get some estimate of which packages need the >> most work. But if there is some area of NumPy or SciPy you are >> familiar with, please go ahead and pitch in. Here is the current >> NumPy/ SciPy coding standard including docstring standards: >> http://projects.scipy.org/scipy/numpy/wiki/CodingStyleGuidelines >> > > I have some questions regarding Travis' latest modifications to the > documentation guidelines: > > The following section was removed, why? > > """ > A reST-documented module should define:: > > __docformat__ = 'restructuredtext en' > > at the top level in accordance with `PEP 258 > <http://www.python.org/dev/peps/pep-0258>`__. Note that the > ``__docformat__`` variable in a package's ``__init__.py`` file does > not apply to objects defined in subpackages and submodules. > """ > I don't see the point in every file having a __docformat__ line, when our documentaion formatting tool should already know the standard we are using is. It's just more cruft. Besides the PEP was rejected, so I don't know why we should be following it. We obviously need a pre-processor to map our files to epydoc, so let's do that instead of contorting docstrings into a format demanded by the tool. The tool is a "tool" for us to use not be chained by. > We had a long discussion on the mailing list on the pros and cons of > "*Parameters*:" vs. "Parameters:". I see now that it has been changed > to > > Parameters > ---------- > > Is this still recognised as a list? > > I noted that the examples are now no longer indented: does ReST allow this? > > > Note that building the example documentation, `epydoc example.py`, now > warns: > > File /tmp/example.py, line 19, in example.foo > Warning: Line 24: Wrong underline character for heading. > Warning: Lines 27, 30, 32, 37, 39, 41, 43, 48, 50: Improper paragraph indentation. > > > While I understand the rationale behind > > "The guiding principle is that human readers of the text itself are > given precedence over contorting the docstring so that epydoc_ > produces nice output." > > I think that it would be impractical to break compatibility with the > only documentation builder we currently have (unless there are others > I am unaware of?). > > In my mind, the rationale trumps the "impracticality" especially since epydoc still produces readable output. It is only the formatting that is not pretty-printed (although it doesn't look bad there either). So, AFAIK, nothing is actually broken (except a little HTML formatting on output of epydoc). But, now the docstrings don't look crufty to my eyes. Basically, the problem is that I *really* don't like the *Parameters*: thing and all the indentation that takes place inside of docstrings to satisfy epydoc. It is also inconsistent with the docstring standard that Enthought uses for the ETS tools suite. It was really hard for me to start writing docstrings that followed that pattern. It seems better to have docstrings than to have them fit into an epydoc_-defined pattern. I'm sorry if I made changes unilaterally. I didn't think there would be much concern as long as docstrings were moving forward. -Travis

2 1

Scikits and stuff
by Travis E. Oliphant 28 Dec '07

28 Dec '07

Hey everyone, In preparation for doc-day tomorrow, I've been thinking about scikits and its relationship to scipy. There seem to be three distinct purposes of scikits: 1) A place for specialized tool boxes that build on numpy and/or scipy that live under a common name-space, but are too specialized to live in scipy itself. 2) A place for GPL or similarly-licensed tools so that people can understand what they are getting and appropriate use can be made. 3) A place for modularly-installed scipy tools. Given these three purposes. It seems that we should have three name-spaces: 1) scikits 2) ??? perhaps scigpl 3) scipy (why not use the same namespace --- there is a technological hurdle that may prevent egg distribution of these until we fix it. But, I think we can fix it and would rather not invent another name-space for something that should be scipy). This idea occurred to me, when I was thinking about "tagging" modules in scikits so users could understand the tool. But, at that point it seemed obvious that we should just have different name-spaces which would promote the needed clarity. Ideas. -Travis O.

10 23