Mailman 3 February 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

rgb_to_hsv in scipy.misc ? (was: [Numpy-discussion] Optimizing speed for large-array inter-element algorithms (specifically, color space conversion))
by Nicolas Pinto 08 Jul '09

08 Jul '09

Hello, Would it be possible to include the following rgb to hsv conversion code in scipy (probably in misc along with misc.imread, etc.) ? What do you think? Thanks in advance. Best regards, -- Nicolas Pinto Ph.D. Candidate, Brain & Computer Sciences Massachusetts Institute of Technology, USA http://web.mit.edu/pinto # ------------------------------------------------------------------------------ import numpy as np def rgb_to_hsv_arr(arr): """ fast rgb_to_hsv using numpy array """ # adapted from Arnar Flatberg # http://www.mail-archive.com/numpy-discussion@scipy.org/msg06147.html # it now handles NaN properly and mimics colorsys.rgb_to_hsv output arr = arr/255. out = np.empty_like(arr) arr_max = arr.max(-1) delta = arr.ptp(-1) s = delta / arr_max s[delta==0] = 0 # red is max idx = (arr[:,:,0] == arr_max) out[idx, 0] = (arr[idx, 1] - arr[idx, 2]) / delta[idx] # green is max idx = (arr[:,:,1] == arr_max) out[idx, 0] = 2. + (arr[idx, 2] - arr[idx, 0] ) / delta[idx] # blue is max idx = (arr[:,:,2] == arr_max) out[idx, 0] = 4. + (arr[idx, 0] - arr[idx, 1] ) / delta[idx] out[:,:,0] = (out[:,:,0]/6.0) % 1.0 out[:,:,1] = s out[:,:,2] = arr_max # rescale back to [0, 255] out *= 255. # remove NaN out[np.isnan(out)] = 0 return out

4 5

Patching fmin_l_bfgs_b in scipy.optimize
by Gilles Rochefort 16 Apr '09

16 Apr '09

Hi, I had some troubles getting a stable version of fmin_l_bfgs_b. I got a segmentation fault as soon as I turned a different value than -1 to iprint argument. As far as I understood, this iprint parameter is directly related to some fortran printing/writing logs. It prints iteration, norm of gradient on screen and also in a file called iterate.dat .. Even recompiling the whole stuff from lapack, atlas etc.. (using gfortran only) up to scipy doesn't solve the problem. Because I definitely gave up the idea to fix the bug by recompiling, I decided to patch the way python interacts with lbfgsb fortran procedure. So, I disabled the fortran iprint and add some iprint functionnality in pure python code directly in lbfgsb.py (scipy/optimize). Doing so, I get rid of the useless iterate.dat file too. While patching, I followed by adding some few things : - an optional stopping rule based on maximum number of iterations. I think this is more usefull than the number of function evaluations because the fortran code performs a linesearch procedure inside the lbfgsb ones. - a callback procedure which allows to inspect solution, criterion, and gradient at current iterate. Maybe my patch should be usefull for someone, so I decided to share it with the scipy community. Gilles. PS: I provide two files : lbfgsb.py -- the full python script to put in scipy/optimize scipy-optimize-lbfgsb-12162008.patch -- a patch file

3 4

Warning about remaining issues in stats.distributions ?
by josef.pktd＠gmail.com 16 Mar '09

16 Mar '09

There are still some known failures or methods, that don't work very well for some distributions: * skew. kurtosis: a few remaining incorrect results, no exact test available * entropy: returns nan in some cases, correctness of values not tested * fit: does not converge, or does not estimate correctly for some distributions, this is expected since maximum likelihood does not always work * distributions that have problems for some range of parameters * incorrect random number generator for one discrete distribution Is it useful to put a warning about these issues in stats.distributions in the release? A friend of mine is doing value-at-risk calculations for operational risk analysis, but I am still reluctant to recommend scipy.stats, when it is not clear whether the results are correct or not, although scripting in Python would be much easier than to struggle with some commercial package with a GUI user interface. The problem is that in this application they are using many distributions automatically, and not just a few that are well tested by frequent use. Per Brodtkorb has an improved version of distributions.py with an additional estimation procedure, which we can discuss as an enhancement after 0.7 is out of the way, if there is an interest in this. Josef

4 15

Scipy workflow (and not tools).
by Matthew Brett 16 Mar '09

16 Mar '09

Hi, I've split this off into a new thread because I felt there were two issues in Stefan's original thread. This is in the hope that we can stimulate discussion on the workflow (as opposed to - say - which version control system to use, or which bugtracker). I would be very interested to see if we can come to a consensus on the important discussion of whether to introduce fairly formal code review into the scipy workflow. I've appended the key piece of discussion below. > 2009/2/24 Travis E. Oliphant <oliphant(a)enthought.com>: >> I think the biggest problem has been time and adding too formal of a >> process will just increase the time it takes to get code into SciPy. >> I'm fine with emphasizing documentation and tests as we discuss things >> and we should encourage each other, but I'm not comfortable with >> hard-line statements like the ones being made above. Yes, such things >> are helpful, but they are also expensive and I worry more about what we >> lose in contributions. > > Having so little time means that we cannot be cavalier about adding > broken code to SciPy. Like Matthew mentioned, this becomes an immense > maintenance burden. > >> The quality of what we create should emerge as all interested parties >> critically look at the code that is available in SciPy. > > I agree with that sentiment; and looking critically at code in SciPy > starts with our own patches. > >> Not everyone >> can do that on the same schedule. I'm opposed to trying to force that >> to happen. I very much favor cultivating a culture that wants someone >> to fix the problems in their code. > > Sure, let's be inclusive, but also set a bar. If you make the time to > write a patch, make the time to do it well (it does not take long to > construct a test -- you have to make sure your code works properly > anyhow). > >> But, my favorite workflow is a bit more chaotic, than that. People >> create their own DVCS versions of SciPy using their best judgment and >> publish revisions they consider to be working code. >> >> Branches that are given the thumbs up by 2 people (or 1 on the steering >> committee) get pushed to the main branch. This review happens >> regularly, on IRC channels at regularly scheduled times. > > Two eyes on every piece of code in SciPy, that's all we need. Two > critical eyes that realise the value of tests and documentation. Your > outline above fits in with my view of how this could happen. Best, Matthew

25 96

talkbox scikit
by Georg Holzmann 04 Mar '09

04 Mar '09

Hallo David ! I want to contribute some code to the talkbox scikit and have some questions. - First, is this scikit also for audio signal processing / audio (music) feature extraction, or mainly for speech only ? - Second, I have implemented some (random) code for audio signal processing which IMHO would be nice to have in a scikit: * Implementation of a Generalized Cross Correlation (GCC) with various pre-whitening filters. (after "The Generalized Correlation Method for Estimation of Time Delay" by Charles Knapp and Clifford Carter, programmed with looking at the matlab GCC implementation by Davide Renzi) this function is used for robustly determine the time delay between two real signals * Equivalent Rectangular Bandwidth Filter Coefficients for biquad IIR Filters. (implemented after "An Efficient Implementation of the Patterson-Holdsworth Auditory Filter Bank" by Malcolm Slaney) * Filter coefficients for a bank of Gammatone filters. (implemented after "An Efficient Implementation of the Patterson-Holdsworth Auditory Filter Bank" by Malcolm Slaney) Implementation also with multiple biquad filters, to avoid numerical unstabilities. * Common filter parameters for audio biquad IIR filters (after "Cookbook formulae for audio EQ biquad filter coefficients", http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt) * Conversion of linear IIR filter parameters to a minimum phase filter with the same amplitude response. * MFCC feature extraction (but I have seen that you already have implemented mfccs...) * I plan to implement more audio/music feature extraction methods in near future (chroma features, beat features, beat-synchronous features ...) - OK, would this be suited for talkbox ? If yes, are there any guidlines how to contribute (I will adapt the code to your scikit style, writing more tests and docs) - I mean, should I just send you the code by mail ? - In which categories should I put all these ? So I propose all the filter parameter calculations in talkbox/fbanks/, feature extraction methods of course into talkbox/features/ and the generalized cross correlation maybe into talkbox/tools/correlations.py, or maybe in a seperate file ... ? - And last but not least, is this the right mailing list for such discussions ;) ? Or are there any special lists for scikits ? OK, thanks for any feedback ! LG Georg

7 14

PEP: Improving the basic statistical functions in Scipy
by Bruce Southey 03 Mar '09

03 Mar '09

Hi, I do apologize in advance if this is considered inappropriate but my goal is to advance the stats capabilities in Scipy. I do recognize that there is a large chunk of excellent work so part of this is simply to ensure that there is adequate documentation and tests. The following is my attempt of a PEP to provide some direction on how to improve the basic statistical functions within Scipy. I do have a list of the individual functions and the arguments involved but I decided it was in appropriate to attach it here. Probably the main aspect that I would like feedback is on whether or not there should be a single interface to these basic statistical functions. Thanks Bruce PEP: Improving the basic statistical functions in Scipy Authors: Bruce Southey Created: 26-Feb-2009 Abstract ======== This current PEP is orientated towards addressing the fundamental problems with the basic statistical functions in Scipy. The outcome is to provide Scipy with a consistent, well-tested and documented set of basic statistical functions that are available to different array types. Motivation ======== This PEP addresses the basic statistical functions available in the stats component of Scipy. These functions are defined in the following files: stats.py – Defines many statistical functions and imports statlib. morestats.py – Adds additional statistical functions to stats.py _support.py - Defines the functions used in stats.py but also is a circular because it also imports stats mstats.py – Just imports functions from mstats_basic.py and mstats_extras.py mstats_basic.py – Defines statistical functions for masked arrays mstats_extras.py– Defines additional statistical functions for masked arrays In total there are 178 unique functions defined in these files, some of which are private or internal and some have the same name but are defined slightly differently between standard and masked arrays. A list of theses functions is available. While the functions are defined for standard arrays and masked arrays, not all functions are available for both array types. For example, the majority of functions defined in stats.py for standard arrays are available in mstats_basic.py. But none of the standard arrays functions defined in morestats.py are available for masked arrays. Also none of these are functions are directly supported for other array types available to Scipy such record arrays, record arrays that contain masked data and sparse arrays. Specification ======== 1) Provide the same basic statistical functions with the same arguments for standard and masked arrays. 2) Utilize a single interface. For example, the gmean function (note _chk_asarray is defined differently): stats.py: def gmean(a, axis=0): a, axis = _chk_asarray(a, axis) log_a = np.log(a) return np.exp(log_a.mean(axis=axis)) mstats_basic.py: def gmean(a, axis=0): a, axis = _chk_asarray(a, axis) log_a = ma.log(a) return ma.exp(log_a.mean(axis=axis)) Rather a single function can be defined as: def gmean(a, axis=0): log_a = np.log(a) return np.exp(log_a.mean(axis=axis)) import numpy as np import numpy.ma as ma X=[1,2,3,4,5] a=np.array(X) m=ma.array(X, mask=[0,0,0,0,0]) np.exp((np.log(X).mean())) #2.6051710846973517 np.exp((np.log(a).mean())) #2.6051710846973517 np.exp((np.log(m).mean())) #2.6051710846973517 3) Depreciation and removal of unnecessary functions such as linregress. 4) Cleanup styles issues including: a) White space usage b) Consistent arguments such as 'a' vs 'x' and the usage of *args c) Uniquely identifying functions. i) rootfunc and tempfunc defined two and three times, respectively, in morestats.py but have different arguments. ii) makestr is defined twice in _support.py, once a main function and once as a subfunction of printcc. 5) Ensure info.py is complete and correct. 6) Improve the documentation of basic statistical functions in connection with the Scipy documentation Marathon (http://www.scipy.org/Developer_Zone/DocMarathon2008) 7) Improve the tests of the basic statistical functions: i) All functions should have at least have basic test coverage that indicates whether or not it is functional. ii) Important functions should have tests that include unexpected elements like Nan's, positive and negative infinity and other unexpected inputs. iii) Ideally there should be tests that check the function accuracy. 8) Extension of the functions to other array types available to Scipy such record arrays, record arrays that contain masked data and sparse arrays. Perhaps beyond the scope of this PEP. Backwards Compatibility ======== There is no guarantee that the outcome will maintain complete backwards compatibility because a consistent API is required across different array types. However, any changes to existing APIs must be justified such as ensuring the same keywords between functions for different array types.

7 24

Improving the bug tracking workflow: starting document
by David Cournapeau 01 Mar '09

01 Mar '09

Hi, Following the discussions, I have started to write a small document highlighting my current gripes with trac. I focus on some common scenario, and pin-point trac limitations. I mention possible new tools at the end, but that's not the main point: everybody who is also disatisfied with trac, and maybe even more importantly people who are currently satisfied and think their scenario is not covered should feel free to comment/modify it: http://scipy.org/scipy/numpy/wiki/ImprovingIssueWorkflow I put the initial version in svn as well: http://projects.scipy.org/scipy/numpy/browser/trunk/doc/neps/newbugtracker.…. cheers, David

5 9

Matrix exponential
by jason-sage＠creativetrax.com 28 Feb '09

28 Feb '09

John Cremona posted the following message to the sage development list about matrix exponentials. I'm copying it to here since it asks about the scipy matrix exponential method (we say numpy below, but we really mean scipy...) John Cremona wrote: > >> I have just been to a colloquium talk by numerical analyst Nick Higham > >> (Manchester) called "How to compute and not to compute a matrix > >> exponential". He has new methods which are now in mathematica, matlab > >> and NAG but (apparantly) nowhere else. He only seemed interested in > >> getting good speed & precision to 16 decimals but (when I asked) > >> confirmed that the methods should apply to give arbitrary precision. > >> > >> I just checked and see that Sage's matrix exp() uses something stupid > >> except over RDF/CDF where it uses a pade approximation method via > >> numpy. The method of the talk was a variant of that, the main trick > >> being to use exactly the right order of Pade approx. so maximise > >> precision and speed. > >> > >> I would like to know how good the numpy method is, and whether it can > >> be improved to this "state of the art" version at least for RDF. Then > >> it could be another selling point for Sage. > > > > Could you CC the numpy devlist as well on this? It sounds exciting! I will if you give me the address (or you can perhaps?). It might be worth including Higham's URL: http://www.maths.manchester.ac.uk/~higham/ as he has lots of his talks up there including some which are similar to the one I heard.

2 3

The future of SciPy and its development infrastructure
by Stéfan van der Walt 28 Feb '09

28 Feb '09

*[If you only have 30 seconds to read this email, read the **bold text only] * *Dear* SciPy *developer*s The past while has seen a rocky ride with the SciPy servers, but yesterday Peter Wang announced that he is attending to the situation. This, then, seems like the perfect time to *stand back and take a look at our infrastructure*, and whether we should continue with the current setup. To put this conversation into context, we have to face the facts: SciPy has a large user community relative to the number of developers. A big library of code, used by many scientists, is supported by a small handful of people all over the world. *We cannot afford* *a high barrier to contribution*, and we have to lower the effort it takes for a developer to merge contributed code. *I'd like to propose two changes* to the status quo: 1. *Change to a distributed revision control system*, encouraging more open collaboration. 2. *Determine guidelines for code acceptance*, in terms of unit tests, documentation and peer review. Allow me to motivate these changes, and then suggest practical approaches for their implementation: Subversion allows only a selected group of developers to change the SciPy source code. This does not encourage a culture of meritocracy, but worse, has practical implications, in that users cannot merge their own patches. I won't discuss the advantages of distributed revision control here, but note that it shifts responsibility from the current core developers to contributers; *that benefits us all!* This ties in with my second point: code review. The current developers have access to SVN because they are experienced programmers with knowledge of SciPy's scientific domains of application. We are unable to employ this scarce resource fully, because it simply takes too long to merge a patch from Trac, review it, *bring it up to scratch*, and commit it. *We have to put a system in place which allows contributers to take responsibility for their own patches, and for core developers to guide and advise during this process.* As it is, we have many patches waiting on Trac for up to a year or more without any feedback; that is not acceptable. My view on testing is simple: *untested code is probably broken code* (and I can show examples from the past year's commit logs to corroborate this statement). *As for documentation, we cannot afford to be without it. * Implementation: Enthought generously hosts SciPy, and I hope they will continue doing so. New software will need to be installed on the server, but we have many hands willing to tackle that task: David Cournapeau and myself included. Before deploying to scipy.org, *we will configure a *different* server as a proof of concept.* 1) *Distributed revision control system: David Cournapeau and myself have been test driving Git [1] on SciPy and NumPy for a while. It is fast, well supported, has great branch support, and is simple to use for the average contributor, while allowing powerful patch-carving for the more adventurous. * 2) *Ticketing back-end:* David is exploring RedMine [2], and I'd like to take a look at InDefero [3], but *we'll do a careful analysis* of trac-git (like FedoraHosted) too. Thank you for taking the time to deliberate on SciPy's future. I would love to hear your comments. Kind regards Stéfan [1] http://git.or.cz/course/svn.html [2] http://www.redmine.org/ [3] http://scipy.indefero.net/p/numpy/ [4] http://fedorahosted.org

36 151