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

prepending underscores to private module names
by Ralf Gommers 16 May '12

16 May '12

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 14 May '12

14 May '12

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) 30 Apr '12

30 Apr '12

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 26 Apr '12

26 Apr '12

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 26 Apr '12

26 Apr '12

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 26 Apr '12

26 Apr '12

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 26 Apr '12

26 Apr '12

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 22 Apr '12

22 Apr '12

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 21 Apr '12

21 Apr '12

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