Mailman 3 July 2009 - 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

Continuous Wavelet Transform for SciPy
by Sean Arms 12 Nov '09

12 Nov '09

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

SciPy Foundation
by Joe Harrington 16 Aug '09

16 Aug '09

About sixteen months ago, I launched the SciPy Documentation Project and its Marathon. Dozens pitched in and now numpy docs are rapidly approaching a professional level. The "pink wave" ("Needs Review" status) is at 56% today! There is consensus among doc writers that much of the rest can be labeled in the "unimportant" category, so we're close to starting the review push (hold your fire, there is a web site mod to be done first). We're also nearing the end of the summer, and it's time to look ahead. The path for docs is clear, but the path for SciPy is not. I think our weakest area right now is organization of the project. There is no consensus-based plan for improvement of the whole toward a stated goal, no centralized coordination of work, and no funded work focused on many of our weaknesses, notwithstanding my doc effort and what Enthought does for code. I define success as popular adoption in preference to commercial packages. I believe in vote-with-your-feet: this goal will not be reached until all aspects of the package and its presentation to the world exceed those of our commercial competition. Scipy is now a grass roots effort, but that takes it only so far. Other projects, such as OpenOffice and Sage, don't follow this model and do produce quality products that compete with commercial offerings, at least on open-source platforms. Before we can even hope for that, we have to do the following: - Docs - Rest of numpy reference pages reviewed and proofed or marked unimportant - Scipy reference pages - User manual for the whole toolstack - Multiple commercial books - Packaging - Personal Package Archive or equivalent for every release of every OS for the full toolstack (There are tools that do this but we don't use them. NSF requires Metronome - http://nmi.cs.wisc.edu/ - for funding most development grants, so right now we're not even on NSF's radar.) - Track record of having the whole toolstack installation "just work" in a few command lines or clicks for *everyone* - Regular, scheduled releases of numpy and scipy - Coordinated releases of numpy, scipy, and stable scikits into PPA system - Public communication - A real marketing plan - Executing on that plan - Web site geared toward multiple audiences, run by experts at that kind of communication - More webinars, conference booths, training, aimed at all levels - Demos, testimonials, topical forums, all showcased - Code - A full design review for numpy 2.0 - No more inconsistencies like median(), lacking "out", degrees option for angle functions? - Trimming of financial functions, maybe others, from numpy? - Package structure review (eliminate "fromnumeric"?) - Goal that this be the last breakage for numpy API (the real 1.0) - Scipy - Is it maintainable? should it be broken up? - Clear code addition path (or decide never to add more) - Docs (see above) - Add-on packages - Both existence of and good indexing/integration/support for field-specific packages - Clearer development path for new packages - Central hosting system for packages (svn, mailing lists, web, build integration, etc.) - Simultaneous releases of stable packages along with numpy/scipy I posted a basic improvement plan some years back. The core ideas have not changed; it is linked from the bottom of http://scipy.org/Developer_Zone. I chose our major weakness to begin with and started the doc project, using some money I could justify spending simply for the utility of docs for my own research. I funded the work of two doc coordinators, one each this summer and last. Looking at http://docs.scipy.org/numpy/stats/, you can see that when a doc coordinator was being paid (summers), work got done. When not, then not. Without publicly announcing what these guys made, I'll be the first to admit that it wasn't a lot. Yet, those small sums bought a huge contribution to numpy through the work of several dozen volunteers and the major contributions of a few. My conclusion is that active and constant coordination is central to motivating volunteer work, and that without a salary we cannot depend on coordination remaining active. On the other hand, I have heard Enthought's leaders bemoan the high cost of devoting employee time to this project, and the low returns available from selling support to universities and non-profit research institutes. Their leadership has moved us forward, particularly in the area of code, but has not provided the momentum necessary to carry us forward on all fronts. It is time for the public and education sectors to kick in some resources and organizational leadership. We are, after all, benefitting immensely. Since the cost of employee time is not so high for us in the public and education sectors, I propose to continue hiring people like Stefan and David as UCF employees or contractors, and to expand to hiring others in areas like packaging and marketing, provided that funding for those hires can be found. However, my grant situation is no longer as rich as it has been the past two years, and the needs going forward are greater than in the past if we're now to tackle all the points above. So, I will not be hiring another doc guru from my research grants next year. I am confident that others are willing to pitch in financially, but few will pitch in a full FTE, and we need several. We can (and will) set up a donations site, but donation sites tend to receive pizza money unless a sugar daddy comes along. Those benefitting most from the software, notably education, non-profit research, and government institutions, are *forbidden* from making donations by the terms of their grants. NSF doesn't give you money so you can give it away. We need to provide services they can buy on subcontract and a means for handling payments from them. Selling support does not solve the problem, as that requires spending most of the income on servicing that particular client. Rather, we need to sell a chunk of documentation or the packaging of a particular release, and then provide the product not just to that client but to everyone. We can also propose directly for federal and corporate grant funds. I have spoken with several NASA and NSF program managers and with Google's Federal Accounts Representative, and the possibilities for funding are good. But, I am not going to do this alone. We need a strong proposal team to be credible. So, I am seeking a group that is willing to work with me to put up the infrastructure of a funded project, to write grant proposals, and to coordinate a financial effort. Members of this group must have a track record of funded grants, business success, foundation support, etc. We might call it the SciPy Foundation. It could be based at UCF, which has a low overhead rate and has infrastructure (like an HR staff), or it might be independent if we can find a good director willing to devote significant time for relatively low pay compared to what they can likely make elsewhere. I would envision hiring permanent coordinators for docs, packaging, and marketing communications. Enthought appears to have code covered by virtue of having hired Travis, Robert, etc.; how to integrate that with this effort is an open question but not a difficult one, I think, as code is our strongest asset at this point. I invite discussion of this approach and the task list above on the scipy-dev(a)scipy.org mailing list. If you are seeing this post elsewhere, please reply only on scipy-dev(a)scipy.org. If you are eligible to lead funding proposals and are interested in participating in grant writing and management activities related to work in our weak areas, please contact me directly. Thanks, --jh-- Prof. Joseph Harrington Planetary Sciences Group Department of Physics MAP 414 4000 Central Florida Blvd. University of Central Florida Orlando, FL 32816-2385 jh(a)physics.ucf.edu planets.ucf.edu

14 38

nbinom.ppf
by Pierre GM 01 Aug '09

01 Aug '09

All, I'm puzzled by this result: >>> from scipy.stats.distributions import nbinom >>> nbinom(.3,.15).ppf(nbinom(.3,.15).cdf(np.arange(20))) array([ 1., 2., 2., 3., 5., 5., 6., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.]) I would have naturally expected np.arange(20). Using np.round instead of np.ceil in nbinom_gen._ppf seems to solve the issue. Is there any reason not to do it ? Thx in advance. P.

3 10

Agenda items for the next steering committee meeting
by Stéfan van der Walt 30 Jul '09

30 Jul '09

Hi all, We are soliciting agenda items for the SciPy Steering Committee meeting taking place during the SciPy09 conference. Please view the related post on the numpy-discussion list and leave any comments there: http://thread.gmane.org/gmane.comp.python.numeric.general/31845 Kind regards Stéfan

1 0

Scipy.org down
by Pauli Virtanen 29 Jul '09

29 Jul '09

Scipy.org website seems to be down. Port 80 is open, but GET yields no response. -- Pauli Virtanen

3 2

Nuther update and no Skypecon this week?
by David Goldsmith 28 Jul '09

28 Jul '09

"Other Math" is now "Goal Met!" Recall (from a previous post to this list) that I can't do the Skypecon tomorrow; no one responded to my request for suggestions of other times, so...I'm assuming no one wants to have one this week (either). DG

1 0

Re: [SciPy-dev] Bug?
by David Goldsmith 28 Jul '09

28 Jul '09

Okee-dokee, thanks! DG --- On Tue, 7/28/09, Ralf Gommers <ralf.gommers(a)googlemail.com> wrote: > From: Ralf Gommers <ralf.gommers(a)googlemail.com> > Subject: Re: [SciPy-dev] Bug? > To: "SciPy Developers List" <scipy-dev(a)scipy.org> > Date: Tuesday, July 28, 2009, 5:58 PM > > > On Tue, Jul 28, 2009 at 8:19 PM, > David Goldsmith <d_l_goldsmith(a)yahoo.com> > wrote: > > > > >>> np.set_printoptions(precision=4) > > >>> np.lib.scimath.arctanh([1j]) > > array([ 0.+0.7854j]) > > > > so far, so good > > > > >>> np.lib.scimath.arctanh(1j) > > 0.78539816339744828j > > > > i.e., set_printoptions didn't affect the result when > the input was a scalar. "Feature" or > bug? > A bug, reported last year in numpy ticket #817. > > Ralf > > > > > > > DG > > > > > > > > _______________________________________________ > > Scipy-dev mailing list > > Scipy-dev(a)scipy.org > > http://mail.scipy.org/mailman/listinfo/scipy-dev > > > > > -----Inline Attachment Follows----- > > _______________________________________________ > Scipy-dev mailing list > Scipy-dev(a)scipy.org > http://mail.scipy.org/mailman/listinfo/scipy-dev >

1 0

Bug?
by David Goldsmith 28 Jul '09

28 Jul '09

>>> np.set_printoptions(precision=4) >>> np.lib.scimath.arctanh([1j]) array([ 0.+0.7854j]) so far, so good >>> np.lib.scimath.arctanh(1j) 0.78539816339744828j i.e., set_printoptions didn't affect the result when the input was a scalar. "Feature" or bug? DG

2 1

Re: [SciPy-dev] random.set_state also in need of EXPERT attention
by David Goldsmith 28 Jul '09

28 Jul '09

Awesome, Pauli, thanks! DG --- On Tue, 7/28/09, Pauli Virtanen <pav+sp(a)iki.fi> wrote: > From: Pauli Virtanen <pav+sp(a)iki.fi> > Subject: Re: [SciPy-dev] random.set_state also in need of EXPERT attention > To: scipy-dev(a)scipy.org > Date: Tuesday, July 28, 2009, 2:44 PM > > (Disclaimer: I'm not a PRNG expert.) > > On 2009-07-28, David Goldsmith <d_l_goldsmith(a)yahoo.com> > wrote: > > Hi, folks! Ralf originally requested expert > attention for this > > in its Discussion section; I over-confidently said I > thought I > > could handle it; after too many hours of unproductive > research > > (I seem to have found that there's a lot of "politics" > in the > > field of random number generators? > > Looks like science as usual to me :) I think the natural > explanation is that everyone proposing a new PRNG of course > wants > to advertise its merits. Also, as there is no clear-cut way > to > quantify the "randomness" of a PRNG and because of speed > tradeoffs, there is something to argue about and many > different > alternatives have been proposed. > > > Despite widespread use: > > http://en.wikipedia.org/wiki/Mersenne_twister, > Matsumoto & > > Nishimura and this algorithm don't appear to have much > > > recognition in certain circles, e.g., those under the > influence > > of Marsaglia, e.g., notably, "Numerical Recipes." > > The second edition of Numerical recipes was written in > 1992, so > that explains why it's not there. > > If you Google it, for the third edition the authors of NR > respond > [1] that they didn't include MT because it "has just too > many > operations per random value generated". This is of course > understandable in a book that aims to give a focused > introduction > on the subject. Whether it reflects the merits of the > algorithm > is then a different question. > > .. [1] http://www.nr.com/forum/showthread.php?t=1724 > > > I understand that Marsaglia is critical of MT, but for > "NR" to > > completely ignore an algorithm in such widespread use, > well, > > out of curiosity, if anyone knows, "what gives"?) > > As I see it, the statements sourced in the Wikipedia > article are > fairly mild, criticising mostly the complexity of the > algorithm. > Also, they were not backed by anything, and probably should > be > taken with a grain of salt. > > Marsaglia seems to have proposed another types of PRNGs in > 2003, > but these had flaws, which maybe were addressed by Brent > later > on. (Cf. [2] and follow the references.) > > The MT article is widely cited (123 citations as reported > by ACM, > 1346 by Google Scholar), and sampling some of the > review-type > ones (eg. [3,4]) the generator seems to have done > reasonably in > various randomness tests and also be reasonable speed-wise. > > Perhaps the algorithm is not optimal -- after all, it's > already > more than ten years old -- but it appears to be well tested > and > understood. > > .. [2] http://wwwmaths.anu.edu.au/~brent/random.html > .. [3] http://dx.doi.org/10.1016/j.csda.2006.05.019 > .. [4] http://www.iro.umontreal.ca/~simardr/testu01/tu01.html > > -- > Pauli Virtanen > > _______________________________________________ > Scipy-dev mailing list > Scipy-dev(a)scipy.org > http://mail.scipy.org/mailman/listinfo/scipy-dev >

1 0