Mailman 3 November 2012 - SciPy-Dev

Bootstrap confidence limits code
by Constantine Evans 05 Oct '14

05 Oct '14

Hello everyone, A few years ago I implemented a scikit for bootstrap confidence limits (https://github.com/cgevans/scikits-bootstrap). I didn’t think much about it after that until recently, when I realized that some people are actually using it, and that there’s apparently been some talk about implementing this functionality in either scipy.stats or statsmodels (I should thank Randal Olson for discussing this and bringing it to my attention). As such I’ve rewritten most of the code, and written up some docstrings. The current code can do confidence intervals with basic percentile interval, bias-corrected accelerated, and approximate bootstrap confidence methods, and can also provide bootstrap and jackknife indexes. Most of it is implemented from the descriptions in Efron and Tibshirani’s Introduction to the Bootstrap, but the ABC code at the moment is a port from the modified-BSD-licensed bootstrap package for R (not the boot package) as I’m not entirely confident in my understanding of the method. And so, I have a few questions for everyone: * Is there any interest in including this sort of code in either scipy.stats or statsmodels? If so, where do people think would be the better place? The code is relatively small; at the moment it is less than 200 lines, with docstrings probably making up 100 of those lines. * Also, if so, what would need to be changed, added, and improved beyond what is mentioned in the Contributing to Scipy part of the reference guide? I’m never a fan of my own code, and imagine quite a bit would need to be fixed; I know tests will need to be added too. In addition, I have a few questions about what would be better practice for the API, and I haven’t really found a guide on best practices for Scipy: * When I started writing the code, I wrote a single function ci for confidence intervals, with a method argument to choose the method. This is easy for users, especially so that they don’t have to look through documentation to realize that BCA is the most generally useful method (at least from everything I’ve read) and that there really isn’t any reason to use many of the simpler methods. However, ABC takes different paramenters, and needs a statistic function that takes weights, which makes this single-function organization trickier. At the moment, I have a separate function for ABC. Would it be better to split up all the methods to their own functions? * ABC requires a statistic function that takes weights. I’ve noticed that things like np.average takes a weights= argument. Would it be better to require input of a stat(data,weights) function, or input of a stat(data,weights=) with weights as a named argument? The latter would be nice in terms of allowing the same function to be used for all methods, but would make it impossible to use a lambda for the function. Is there some other method of doing this entirely? * Are there any missing features that anyone thinks should be added? I apologize if much of this is answered elsewhere, I just haven’t found any of it; I also apologize if this is far too long-winded and confusing! Regards, Constantine Evans

4 3

SLSQP Constrained Optimizer Status
by Rob Falck 25 Nov '13

25 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

Scipy.org is down + question about "topical software"
by Michael Clerx 03 Jan '13

03 Jan '13

Dear all, First of all, scipy.org appears to be down and has been for some time today. Yesterday, when it was still up, I registered an account and tried adding a link to the topical software section but got the error "too many links". Does anybody know how/if I can fix this? kind regards, Michael

7 12

PR: Add solout support to dopri5 and dop853
by John Travers 06 Dec '12

06 Dec '12

Hi all, I realise that the ode code has its problems/detractors, but I needed a long requested feature, so I tried to add it: https://github.com/scipy/scipy/pull/366 This enables the use of a callback function which is called at each internal integration step of dopri5 or dop853. This is useful to monitor the ode evolution at its natural time-scale without increasing the output density. It also enables stopping the integration mid-way. This feature was requested quite a while ago, just before my absence from scipy development: http://mail.scipy.org/pipermail/scipy-dev/2010-September/015540.html Comments/improvements welcome. Cheers, John

2 2

Gmres with Sparse Complex Matrices gives Segmentation fault (or NULL result)
by Martin Campos Pinto 04 Dec '12

04 Dec '12

Hi all, I am having some troubles using iterative solvers with sparse complex matrices. (I first posted that issue on the SciPy-User mailing list but I guess this one is more appropriate.) So, here is an example with gmres: the script from scipy import sparse from numpy.linalg import norm from numpy.random import rand import scipy.sparse.linalg as spla C = sparse.lil_matrix((10, 10), dtype=complex) C.setdiag(rand(10)) C[0,3] = 1j C = C.tocsr() c = rand(10)+rand(10)*1j x = spla.spsolve(C, c) print "spsolve: norm(C*x - c) = ", norm(C*x - c) (x,info) = spla.gmres(C, c) print "gmres: norm(C*x - c) = ", norm(C*x - c) gives as output: spsolve: norm(C*x - c) = 1.57621370006e-16 Segmentation fault Actually, sometimes I get this error message instead of a Segmentation fault: Traceback (most recent call last): File "test_gmres_cplx.py", line 45, in <module> (x,info) = spla.gmres(C, c) File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/iterative.py", line 437, in gmres revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob) SystemError: NULL result without error in PyObject_Call Note that I have tested gmres with real sparse matrices, and it runs fine: Indeed C = sparse.lil_matrix((10, 10)) C.setdiag(rand(10)) C = C.tocsr() c = rand(10) x = spla.spsolve(C, c) print "spsolve: norm(C*x - c) = ", norm(C*x - c) (x,info) = spla.gmres(C, c) print "gmres: norm(C*x - c) = ", norm(C*x - c) gives spsolve: norm(C*x - c) = 6.93889390391e-18 gmres: norm(C*x - c) = 5.28860261481e-16 I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer, and I imagine that if those were due to bugs, they should have been fixed by now... Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ? (I am using scipy version 0.10.1) Thanks, Martin -- Martin Campos Pinto LJLL, UPMC & CNRS

2 5

Re: [SciPy-Dev] [SciPy-User] Request for Feedback/Authorization: a new public repository for the conference
by Jonathan Rocher 29 Nov '12

29 Nov '12

Thanks to all. Jonathan On Wed, Nov 28, 2012 at 1:17 PM, Pauli Virtanen <pav(a)iki.fi> wrote: > 28.11.2012 20:37, Robert Kern kirjoitti: > > On Wed, Nov 28, 2012 at 7:30 PM, Pauli Virtanen <pav(a)iki.fi> wrote: > >> 28.11.2012 20:28, Robert Kern kirjoitti: > >> [clip] > >>> +1. If there is no strenuous objection in short order, I will make a > >>> "SciPy Conference" team under the scipy organization and add you two > >>> to it. > >> > >> +1 also from me. > > > > On the basis of exit polls, I'm calling this election early. Done. > > The repo is now also there, so things should be ready to roll. > > -- > Pauli Virtanen > > > _______________________________________________ > SciPy-User mailing list > SciPy-User(a)scipy.org > http://mail.scipy.org/mailman/listinfo/scipy-user > -- Jonathan Rocher, PhD Scientific software developer Enthought, Inc. jrocher(a)enthought.com 1-512-536-1057 http://www.enthought.com

1 0

Request for Feedback/Authorization: a new public repository for the conference
by Jonathan Rocher 29 Nov '12

29 Nov '12

Dear all, Andy Terrel and I are the co-chairs of the Scipy2013 conference in Austin next June in Austin. We are very excited about this coming edition and will post regularly to these lists for feedback and announcements. But for now, we would like to consolidate its organization under a public open source project inside the Scipy organization on github. We propose that the project be called scipy-conference. The repo would host documents useful to run the conference as well as code to set up the website and the backend. I write here to ask if the community is ok with that and to ask for the rights to create/manage that new repository. Thanks in advance, Jonathan -- Jonathan Rocher, PhD Scientific software developer Enthought, Inc. jrocher(a)enthought.com 1-512-536-1057 http://www.enthought.com

4 5

Python Scipy Block Diagram Editor with ACG
by dwight reid 26 Nov '12

26 Nov '12

Hey Guys, I am Dwight Reid an electrical engineer and lecturer at the University of Technology (also IEEE member and branch counselor) in Jamaica. Some months ago I had started working on a block diagram editor for python that I dream one day will become a reasonable alternative to Matlab Simulink. I want to put a special emphasis on automatic code generation (ACG) and so far I have done ACG for python Scipy and also for the arduino, these are however just an implementation of some basic functions. As is, the block diagram editing is done in libre office draw and you have to manually type in the functions and commands, I have attached an example exported to pdf. Would anyone be willing to help with this, I have no experience developing with the open source community. I think it would definitely get more people using python since it is much easier to program graphically and it would be very useful to scientific computing. I imagine a better, moretailored graphic editor with predefined blocks that you can pull out like simulink or scilab xcos (maybe a stripped down libre office draw), implementing this is beyond my current skill level so it goes slowly. I hope I'm clear enough, let me know what you think. Dwight

9 14

ANN: SfePy 2012.4
by Robert Cimrman 22 Nov '12

22 Nov '12

I am pleased to announce release 2012.4 of SfePy. Description ----------- SfePy (simple finite elements in Python) is a software for solving systems of coupled partial differential equations by the finite element method. The code is based on NumPy and SciPy packages. It is distributed under the new BSD license. Home page: http://sfepy.org Downloads, mailing list, wiki: http://code.google.com/p/sfepy/ Git (source) repository, issue tracker: http://github.com/sfepy Highlights of this release -------------------------- - initial support for hierarchical basis on quadrilateral and brick elements - unified C/Cython structures for reference mappings - new linear combination boundary condition: edge direction - new examples showing some advanced features For full release notes see http://docs.sfepy.org/doc/release_notes.html#id1 (rather long and technical). Best regards, Robert Cimrman and Contributors (*) (*) Contributors to this release (alphabetical order): Bjarke Dalslet, Vladimír Lukeš, Matyáš Novák

1 0

Use OpenBLAS for the binary releases?
by Sturla Molden 20 Nov '12

20 Nov '12

I think NumPy and SciPy should consider to use OpenBLAS (a fork of GotoBLAS2) instead of ATLAS or f2c'd Netlib BLAS for the binary releases. Here are its virtues: * Very easy to build: Just a makefile, no configuration script or special build tools. * Building ATLAS can be a PITA. So why bother? * Faster than ATLAS, sometimes faster than MKL. * Multithreaded BLAS kernels: OpenMP on Unix, Windows threads on Windows. * The quality of its ancestor GotoBLAS is undisputed. I was the BLAS implementation of choice for major HPC projects around the World. * Free as in BSD licensed. * Funded and developed for use in major Chinese HPC projects. Actively maintained. (GotoBLAS2 is abandonware.) * Open source. The C sources are a pleasure to read, and very easy to verify. * No OpenMP on Windows means no dependency on pthreads-win32 (an LGPL library) when building with MinGW. * Builds on Windows with MinGW and MSYS, and perhaps even without MSYS. * Cygwin is not needed on Windows (this is just BS from the GotoBLAS documentation). Thus, 64-buit builds are possible (I've built it using TDM-GCC for Win64 and 32-bit MSYS). Sturla

1 0