Mailman 3 April 2012 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck Nov. 24, 2013

Nov. 24, 2013

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

prepending underscores to private module names
by Ralf Gommers May 16, 2012

May 16, 2012

Hi, A while ago we discussed making a clear distinction between public and private modules by using underscores consistently, see http://thread.gmane.org/gmane.comp.python.scientific.devel/14936. I've done this for a few modules now, see https://github.com/rgommers/scipy/tree/underscores. Most files could simply be underscored, some other would conflict with extension modules so I also had to append _py to the name. Also I added a file doc/source/api.rst that states which modules are public (up for discussion of course). With google code search I checked how common it is to import from modules that are supposed to be private, for example "from scipy.optimize.optimize import fmin". It's not very common but does happen, so deprecation warnings should be put in. Before going further I'd like to check that this looks okay to people. What do you think? Cheers, Ralf <https://github.com/rgommers/scipy/tree/underscores>

4 4

scipy.stats: algorithm to for ticket 1493
by nicky van foreest May 14, 2012

May 14, 2012

Hi, I saw Ralf's suggestions on issues of scipy.stats to work on. In another project I had to solve a similar problem as the one mentioned in http://projects.scipy.org/scipy/ticket/1493. I like to propose the algorithm below to tackle this. I am looking forward to your feedback. Once the algo seems to work impeccably, I'll try to convert it such that it can be incorporated in scipy.stats and add test cases. Lets first test the example of the ticket: In [21]: from scipy.stats import invnorm In [22]: from scipy import optimize In [23]: invnorm.xb = 10 In [24]: sol = invnorm.ppf(0.8455, 7.24000019602, scale= 2.51913630166) In [25]: print sol 25.1878345282 In [26]: print invnorm.cdf(sol, 7.24000019602, scale=2.51913630166) 0.8455 In [27]: Weird that this leads to a solution, as the solution is not in the search interval. Trying xb = 5 leads to the reported error. Now a solution: # Proposed solution. # A quintessential property in the algorithm is that the cdf is an # non-decreasing function. def findppf(q): # search until xb is large enough left, right = invnorm.xa, invnorm.xb while invnorm.cdf(right, 7.24000019602, scale=2.51913630166) < q: left = right right *= 2 return optimize.brentq(lambda x: \ invnorm.cdf(x, 7.24000019602, scale=2.51913630166) - q,\ left, right) # Perhaps increasing with a larger factor than 2 is faster, as brentq # convergest very fast. Taking `right` too large is problably not a real # problem. On the other hand, when `right` increases too fast, cdf(right) may # become numerically equal to 1. So, what would be a good factor? sol = findppf(0.8455) print sol print invnorm.cdf(sol, 7.24000019602, scale=2.51913630166) This code (I ran it in a file, not in ipython) gives the right results, at least for these values. Nicky

2 21

SciPy Dev Wiki bug Trac broken - how to report bug?
by Pundurs Mark (Nokia-LC/Chicago) April 30, 2012

April 30, 2012

http://projects.scipy.org/scipy/newticket returns an error - with only the unhelpful suggestion to open a ticket: ************** Oops... Trac detected an internal error: OperationalError: database is locked There was an internal error in Trac. It is recommended that you inform your local Trac administrator and give him all the information he needs to reproduce the issue. To that end, you could Create a ticket. The action that triggered the error was: POST: /newticket ************** Is there any other means to register the following bug report? scipy.cluster.hierarchy.ClusterNode.pre_order returns IndexError for non-root node 0.9.0 To reproduce the error, run the following script: import random, numpy from scipy.cluster.hierarchy import linkage, to_tree datalist = [] for i in range(8000): datalist.append(random.random()) datalist = numpy.array(datalist) datalist = numpy.reshape(datalist, (datalist.shape[0], 1)) Z = linkage(datalist) root_node_ref = to_tree(Z) left_root_node_ref = root_node_ref.left left_root_node_ref.pre_order() The result is: Traceback (most recent call last): File "C:\ReproduceError-pre_order.py", line 12, in <module> left_root_node_ref.pre_order() File "C:\Python27\lib\site-packages\scipy\cluster\hierarchy.py", line 732, in pre_order if not lvisited[ndid]: IndexError: index out of bounds One possible solution (successfully tested with preceding script) is to change pre_order in hierarchy.py as follows: n = self.count curNode = [None] * (2 * n) #following two lines changed: dictionaries instead of lists lvisited = {} rvisited = {} curNode[0] = self k = 0 preorder = [] while k >= 0: nd = curNode[k] ndid = nd.id if nd.is_leaf(): preorder.append(func(nd)) k = k - 1 else: #following line changed: check existence of dictionary key rather than value of list item if ndid not in lvisited.keys(): curNode[k + 1] = nd.left lvisited[ndid] = True k = k + 1 #following line changed: check existence of dictionary key rather than value of list item elif ndid not in rvisited.keys(): curNode[k + 1] = nd.right rvisited[ndid] = True k = k + 1 else: k = k - 1 return preorder Mark Pundurs Data Analyst - Traffic Location & Commerce Chicago The information contained in this communication may be CONFIDENTIAL and is intended only for the use of the recipient(s) named above. If you are not the intended recipient, you are hereby notified that any dissemination, distribution, or copying of this communication, or any of its contents, is strictly prohibited. If you have received this communication in error, please notify the sender and delete/destroy the original message and any copy of it from your computer or paper files.

2 1

scipy.stats: some questions/points about distributions.py + reply on ticket 1493
by nicky van foreest April 26, 2012

April 26, 2012

Hi Josef, Sorry for not responding earlier... too many obligations. Before I get back to your earlier mail, I have some naive questions about distributions.py. I hope you don't mind that I fire them at you. 1: https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L436 Is this code "dead"? Within distributions.py it is not called. Nearly the same code is written here: https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L11… 2: I have a similar point about: https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L358 What is the use of this code? It is not called anywhere. Besides this, from our discussion about ticket 1493, this function returns the centralized moments, while the "real" moment E(X^n) should be returned. Hence, the code is also not correct, i.e., not in line with the documentation. 3: Suppose we would turn xa and xb into private atrributes _xa and _xb, then i suppose that https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L883 requires updating. 4: I have a hard time understanding the working (and goal) of https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L883 Where is the right place to ask for some clarification? Or should I just think harder? 5: The definition of arr in https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L60 does not add much (although it saves some characters at some points of the code), but makes it harder to read the code for novices like me. (I spent some time searching for a numpy function called arr, only to find out later that it was just a shorthand only used in the distribution.py module). Would it be a problem to replace such code by the proper numpy function? 6: https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L538 contains a typo. It should be Weisstein. 7: https://github.com/scipy/scipy/blob/master/scipy/stats/distributions.py#L625 This code gives me even a harder time than _argsreduce. I have to admit that I simply don't know what this code is trying to prevent/check/repair. Would you mind giving a hint? Nicky

2 4

Cython 0.16
by Warren Weckesser April 26, 2012

April 26, 2012

Is there currently a constraint on the version of cython that can be used in scipy? If so, can we bump it up to the latest version (0.16)? I would like to take advantage of the fused types in a extension module that I'm working on. Warren

5 5

SciPy 2012 - The Eleventh Annual Conference on Scientific Computing with Python
by Warren Weckesser April 26, 2012

April 26, 2012

SciPy 2012, the eleventh annual Conference on Scientific Computing with Python, will be held July 16–21, 2012, in Austin, Texas. At this conference, novel scientific applications and libraries related to data acquisition, analysis, dissemination and visualization using Python are presented. Attended by leading figures from both academia and industry, it is an excellent opportunity to experience the cutting edge of scientific software development. The conference is preceded by two days of tutorials, during which community experts provide training on several scientific Python packages. Following the main conference will be two days of coding sprints. We invite you to give a talk or present a poster at SciPy 2012. The list of topics that are appropriate for the conference includes (but is not limited to): - new Python libraries for science and engineering; - applications of Python in solving scientific or computational problems; - high performance, parallel and GPU computing with Python; - use of Python in science education. Specialized Tracks Two specialized tracks run in parallel to the main conference: - High Performance Computing with Python Whether your algorithm is distributed, threaded, memory intensive or latency bound, Python is making headway into the problem. We are looking for performance driven designs and applications in Python. Candidates include the use of Python within a parallel application, new architectures, and ways of making traditional applications execute more efficiently. - Visualization They say a picture is worth a thousand words--we’re interested in both! Python provides numerous visualization tools that allow scientists to show off their work, and we want to know about any new tools and techniques out there. Come show off your latest graphics, whether it’s an old library with a slick new feature, a new library out to challenge the status quo, or simply a beautiful result. Domain-specific Mini-symposia Mini-symposia on the following topics are also being organized: - Computational bioinformatics - Meteorology and climatology - Astronomy and astrophysics - Geophysics Talks, papers and posters We invite you to take part by submitting a talk or poster abstract. Instructions are on the conference website: <http://conference.scipy.org/scipy2012/papers.php> http://conference.scipy.org/scipy2012/talks.php<http://conference.scipy.org/scipy2012/papers.php> <http://conference.scipy.org/scipy2012/papers.php> Selected talks are included as papers in the peer-reviewed conference proceedings, to be published online. Tutorials Tutorials will be given July 16–17. We invite instructors to submit proposals for half-day tutorials on topics relevant to scientific computing with Python. See http://conference.scipy.org/scipy2012/tutorials.php<http://conference.scipy.org/scipy2011/tutorials.php> <http://conference.scipy.org/scipy2011/tutorials.php> for information about submitting a tutorial proposal. To encourage tutorials of the highest quality, the instructor (or team of instructors) is given a $1,000 stipend for each half day tutorial. Student/Community Scholarships We anticipate providing funding for students and for active members of the SciPy community who otherwise might not be able to attend the conference. See http://conference.scipy.org/scipy2012/student.php<http://conference.scipy.org/scipy2011/student.php> <http://conference.scipy.org/scipy2011/student.php> for scholarship application guidelines. Be a Sponsor The SciPy conference could not run without the generous support of the institutions and corporations who share our enthusiasm for Python as a tool for science. Please consider sponsoring SciPy 2012. For more information, see http://conference.scipy.org/scipy2012/sponsor/index.php Important dates: Monday, April 30: Talk abstracts and tutorial proposals due. Monday, May 7: Accepted tutorials announced. Monday, May 13: Accepted talks announced. Monday, June 18: Early registration ends. (Price increases after this date.) Sunday, July 8: Online registration ends. Monday-Tuesday, July 16 - 17: Tutorials Wednesday-Thursday, July 18 - July 19: Conference Friday-Saturday, July 20 - July 21: Sprints We look forward to seeing you all in Austin this year! The SciPy 2012 Team http://conference.scipy.org/scipy2012/organizers.php

1 1

scipy.sparse and OpenMP
by Maximilian Nickel April 26, 2012

April 26, 2012

Hi everyone, I've been working with fairly large sparse matrices on a multiprocessor system lately and noticed that scipy.sparse is single-threaded. Since I needed faster computations, I've quickly added some OpenMP #pragma directives in scipy/sparse/sparsetools to the functions that I've been using in order to enable multithreading, what worked out nicely. I wondered if you would be interested in a more complete OpenMP-enabled version of scipy.sparse.setuptools. I've attached the patch of the quick-and-dirty changes that I made so far to this mail, to give you an idea. Best regards Max

5 4

scipy.stats
by nicky van foreest April 22, 2012

April 22, 2012

Hi, While reading the code and the examples of scipy.stats I came across the attribute rv.dist. Specifically, line 227 of the scipy.stats source on the gibhub server mentions: >>> x = np.linspace(0, np.minimum(rv.dist.b, 3)) The problem is that the meaning of this attribute is nowhere explained. Is this a bug? Nicky

2 16

Contributing to SciPy guide
by Ralf Gommers April 21, 2012

April 21, 2012

Hi all, It has been pointed out by a number of people that it's not so easy to get started with contributing to SciPy, and better documentation may help here. So I've written a guide for this. It would be great to get some feedback especially from people who've found it difficult to find this information before. And if you haven't contributed before but were thinking about doing so, perhaps this is a good opportunity to get started! Pull Request: https://github.com/scipy/scipy/pull/191 Rendered guide: https://github.com/rgommers/scipy/blob/howto-contribute/doc/HOWTO_CONTRIBUT… Cheers, Ralf

4 9

My Intro
by Patrick "Kai" Baker April 21, 2012

April 21, 2012

I have never worked on an open source project before so I'm not exactly sure what I'm doing. But I'm real excited about doing this. I have an advanced degree in mathematics and an emphasis on mathematical physics. I love numerical analysis and coding so I'm real excited about joining the SciPy team. At some point, once I get a handle on things, I will like to contribute to the documentation as well. I guess I will hang out, figure out what's going on, read the archives, and check out the source code. Real excited about this! Regards, Kai

4 8

PR #196: simplification of optimization wrappers
by Denis Laxalde April 19, 2012

April 19, 2012

Hi, In pull request #196, I am proposing to simplify the signature of optimization wrappers (minimize, minimize_scalar, etc.) to get something like: x, info = minimize(fun, x0, [jac, constraints], options) This boils down to eliminating the full_output and retall parameters from respective function. Besides, the info dictionnary would always be returned. Comments welcome - https://github.com/scipy/scipy/pull/196 Thanks, -- Denis

6 11

PR #195: optimize.root - interface to root finding algorithms
by Denis Laxalde April 18, 2012

April 18, 2012

Hi, In pull request #195, I am proposing the addition of a unified interface `optimize.root` for root finding algorithms for multivariate functions, similar to `optimize.minimize` and `optimize.minimize_scalar` (already merged). Available solvers are: - hybr: MINPACK's hybrid Powell algorithm (fsolve) - lm: MINPACK's Levenberg-Marquardt algorithm (leastsq) - and solvers from optimize.nonlin (broyden1, krylov, etc.) See <https://github.com/scipy/scipy/pull/195>. Comments welcome. Thanks, Denis.

1 0

testbench for some non-derivative optimizers in NLopt and scipy.optimize
by denis April 17, 2012

April 17, 2012

Folks, fwiw, https://github.com/denis-bz/opt has testbenches for some of the non-derivative optimizers in NLopt and scipy.optimize. >From opt-py.Readme -- These few test functions are *very* noisy, sensitive to random seed: fmin 0.4 neval 369 LN_BOBYQA fmin 3 neval 580 LN_BOBYQA fmin 0.06 neval 423 LN_BOBYQA fmin 12 neval 385 LN_BOBYQA fmin 2 neval 325 LN_BOBYQA ... ==> we need test functions. Comments, any birds of a feather ? cheers -- denis

2 1

ANN: SpikeSort 0.12
by Bartosz Telenczuk April 16, 2012

April 16, 2012

We are pleased to announce the first official release of the SpikeSort, a new spike sorting library based on dynamic and interactive language (Python). SpikeSort 0.12 is available for download at http://spikesort.org SpikeSort aims to be both flexible and user-friendly. It achieves that by defining a set of components that can be mix-and-matched to fit specific needs. These components are based on standard scientific libraries, so they can be easily re-used in custom applications. Main features: * user-friendly and customisable, * interactive command-line interface in Python, * visualization widgets, * k-means and gaussian-mixture-models (GMM) clustering algorithms and manual cluster cutting, * support for multi-channel data (for example, from tetrodes), * support for binary datasets and HDF5 files (support for other formats planned in future version). SpikeSort includes >60 tests making it one of the best-tested spike sorting software on the market. The distribution also contains full documentation, sample scripts and sample data, so that you can start playing with it almost immediately. SpikeSort is offered free of charge under a liberal Open Source license (two-clause BSD license) allowing for non-commercial and commercial use. The project was partially supported by Deutsche Forschungsgemeischaft (SFB 618 ``Theoretische Biologie'', project B4). The sample data is available by courtesy of Stuart Baker and Wellcome Trust who provided funding for data collection. Yours, Bartosz Telenczuk Dmytro Bielievtsov

1 0

Trac unlock?
by Ralf Gommers April 15, 2012

April 15, 2012

Hi, SciPy Trac now seems permanently locked. Can someone who knows how please kick it so it works again after a few tries (or even on the first try)? Thanks, Ralf

2 2

Could anyone compile win32 Py-2.7.2 module for SciPy-0.11dev as is?
by Евгений Черкашин April 12, 2012

April 12, 2012

Hi! I need some last features of signal processing developed in SciPy-0.11dev for merging in my project for Windows users, but I am not able to compile it myself for Win32 and Python-2.7.2. Can Anyone help me and send me a compiled module. It not necessary to the module to be stable one. Thank You and Best regards, Evgeny

1 0

wrapped FITPACK's `sphere.f`
by Andreas H. April 11, 2012

April 11, 2012

Hi, I just wrapped FITPACK's `sphere.f` for scipy.interpolate, to make interpolation/smoothing in spherical coordinates possible. Here's the pull request: https://github.com/scipy/scipy/pull/192 Sorry for not discussing this here earlier, I just now found out about the "Contributing to SciPy" guide. Cheers, Andreas.

1 0

stats.ttest_ind doc
by Deniz Turgut April 6, 2012

April 6, 2012

docs [1] say that ttest_ind: Calculates the T-test for the means of TWO INDEPENDENT samples of scores. This is a two-sided test for the null hypothesis that 2 independent samples have identical average (expected) values. This function performs t-test for two independent samples with *identical variances*. There is also a form of t-test for independent samples with different variances, also known as Welch's t-test [2]. I think it is better to include the 'identical variance' assumption in the doc to avoid confusion. PS: I could do the edit. I registered but looks like I need edit rights. And apparently this is the place to ask for it. My username: avaris [1] http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.h… [2] http://en.wikipedia.org/wiki/Welch's_t_test

2 1

cephes_smirnov never returns on mips/sparc/...
by Yaroslav Halchenko April 1, 2012

April 1, 2012

I have reported this issue some time ago on Debian http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=653948 and then forgot about it until now that I ran into it again. With recent scipy (including git master 0fbfdbc) scipy.stats.ksone.fit seems to stall (never return) on big-endian boxes and return (1.0, nan, nan) on x86. It seems that it is working correctly with scipy 0.7.2 (as it is now in Debian stable). So -- is that anything known? snippet to replicate: import numpy as np import scipy.stats as ss d = np.array([-0.18879233, 0.15734249, 0.18695107, 0.27908787, -0.248649, -0.2171497 , 0.12233512, 0.15126419, 0.03119282, 0.4365294 , 0.08930393, -0.23509903, 0.28231224, -0.09974875, -0.25196048, 0.11102028, 0.1427649 , 0.10176452, 0.18754054, 0.25826724, 0.05988819, 0.0531668 , 0.21906056, 0.32106729, 0.2117662 , 0.10886442, 0.09375789, 0.24583286, -0.22968366, -0.07842391, -0.31195432, -0.21271196, 0.1114243 , -0.13293002, 0.01331725, -0.04330977, -0.09485776, -0.28434547, 0.22245721, -0.18518199, -0.10943985, -0.35243174, 0.06897665, -0.03553363, -0.0701746 , -0.06037974, 0.37670779, -0.21684405]) print "Fitting now" print ss.ksone.fit(d) -- =------------------------------------------------------------------= Keep in touch www.onerussian.com Yaroslav Halchenko www.ohloh.net/accounts/yarikoptic

3 21