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

scipy.test() segfault
by Nils Wagner 06 Nov '12

06 Nov '12

Program received signal SIGSEGV, Segmentation fault. 0x00007ffff7ab2b26 in PyObject_RichCompare () from /usr/lib64/libpython2.7.so.1.0 (gdb) bt #0 0x00007ffff7ab2b26 in PyObject_RichCompare () from /usr/lib64/libpython2.7.so.1.0 #1 0x00007ffff7ab2def in PyObject_RichCompareBool () from /usr/lib64/libpython2.7.so.1.0 #2 0x00007fffe3ba8584 in OBJECT_compare (ip1=<optimized out>, ip2=<optimized out>) at scipy/signal/sigtoolsmodule.c:827 #3 0x00007ffff74b8b64 in msort_with_tmp.part.0 () from /lib64/libc.so.6 #4 0x00007ffff74b8e9b in qsort_r () from /lib64/libc.so.6 #5 0x00007fffe3ba90d6 in PyArray_OrderFilterND (op1=<optimized out>, op2=<optimized out>, order=54989472) at scipy/signal/sigtoolsmodule.c:969 #6 0x00007fffe3ba92a7 in sigtools_order_filterND (__NPY_UNUSED_TAGGEDdummy=<optimized out>, args=<optimized out>) at scipy/signal/sigtoolsmodule.c:1135 #7 0x00007ffff7affcaa in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #8 0x00007ffff7b059c3 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 #9 0x00007ffff7aff7ce in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #10 0x00007ffff7b0204f in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #11 0x00007ffff7b05727 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 #12 0x00007ffff7a9ddc8 in ?? () from /usr/lib64/libpython2.7.so.1.0 #13 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #14 0x00007ffff7b00cae in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #15 0x00007ffff7b05727 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 #16 0x00007ffff7a9dba1 in ?? () from /usr/lib64/libpython2.7.so.1.0 #17 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #18 0x00007ffff7a87c8e in ?? () from /usr/lib64/libpython2.7.so.1.0 #19 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #20 0x00007ffff7acb84a in ?? () from /usr/lib64/libpython2.7.so.1.0 #21 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #22 0x00007ffff7affe75 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #23 0x00007ffff7b0204f in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #24 0x00007ffff7b05727 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 #25 0x00007ffff7a9ddc8 in ?? () from /usr/lib64/libpython2.7.so.1.0 #26 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #27 0x00007ffff7b00cae in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #28 0x00007ffff7b05727 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 #29 0x00007ffff7a9dba1 in ?? () from /usr/lib64/libpython2.7.so.1.0 #30 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #31 0x00007ffff7a87c8e in ?? () from /usr/lib64/libpython2.7.so.1.0 #32 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #33 0x00007ffff7acb84a in ?? () from /usr/lib64/libpython2.7.so.1.0 #34 0x00007ffff7a7d464 in PyObject_Call () from /usr/lib64/libpython2.7.so.1.0 #35 0x00007ffff7affe75 in PyEval_EvalFrameEx () from /usr/lib64/libpython2.7.so.1.0 #36 0x00007ffff7b05727 in PyEval_EvalCodeEx () from /usr/lib64/libpython2.7.so.1.0 >> scipy.__version__ '0.12.0.dev-bb436fa' x86_64 GNU/Linux

3 8

Argument casts in special functions
by Pauli Virtanen 20 Oct '12

20 Oct '12

Hi, Related to: https://github.com/scipy/scipy/pull/334 Currently we have the following behavior for the integer-order Bessel function: from scipy.special import yn print yn(1, 0.123) # -> -5.2816754510498471 print yn(1.5, 0.123) # -> -5.2816754510498471 That is, the floating point argument is silently truncated to an integer. This is not so optimal: the more strict behavior would be to return a `nan` to signify the domain error. The question now is: is it worth preserving the old behavior of the existing functions? My guess would be yes, but I think the unsafe casts should be used just as a legacy special case, and not as a general rule. (Also, the integer-order functions could be deprecated --- there is no speed advantage in them as the kernel selection can be done based on input / input data type). -- Pauli Virtanen

2 1

optimize: add algorithm for global optimization basin hopping
by Jacob Stevenson 18 Oct '12

18 Oct '12

Hi everyone, I just submitted a pull request to incorporate a basinhopping routine in the optimize package. Basinhopping is a powerful algorithm, but all the hard work would be done by the minimizers that already exist in scipy.optimize, so the additional code needed is really not very much. The following is from the documentation notes: Basin hopping is a random algorithm which attempts to find the global minimum of a smooth scalar function of one or more variables. The algorithm was originally described by David Wales http://www-wales.ch.cam.ac.uk/ . The algorithm is iterative with each iteration composed of the following steps 1) random displacement of the coordinates 2) local minimization 3) accept or reject the new coordinates based on the minimized function value. This global minimization method has been shown to be extremely efficient on a wide variety of problems in physics and chemistry. It is especially efficient when the function has many minima separated by large barriers. See the Cambridge Cluster Database http://www-wales.ch.cam.ac.uk/CCD.html for database of molecular systems that have been optimized primarily using basin hopping. This database includes minimization problems exceeding 300 degrees of freedom. Thanks, Jake p.s. this is my first post and first submission

7 8

Making lobpcg support complex matrix
by Gregory Crosswhite 18 Oct '12

18 Oct '12

Dear all, At the moment lobpcg seems to have been designed with only real numbers in mind; this is unfortunate because I have been having some trouble with eigs and was hoping to try it out as an alternative. Fortunately, it looks to me like the solution is rather simple. I was able to get it to work by replacing all instances of ".T" with ".T.conj()" --- which meant that all of the dot products were now correct for complex numbers --- and by replacing the cast to "float64" with a cast to "complex128". Changing the second cast is of course less than ideal in the cases where "float64" is indeed what is being used, but something like it is needed to make the answers not be garbage for complex input matrices. With these changes, I generated 10x10 complex Hermitian matrices A and B to serve as respectively the problem matrix and the normalization/metric matrix, and ran lobpcg with a tolerance of 1e-10. For the standard eigenvalue problem with k=2 (and random X) lobpcg got essentially the same answers as eigh for the lowest two eigenvalue/eigenvector pairs, and for the generalized eigenvalue problem it got the same eigenvalues as eigh and it got eigenvectors whose entries were within ~10^-4 of the eigenvectors returned by eigh (after dividing by the first entry to make the vectors proportionately the same) and such that the normalized overlap between the eigenvectors of lobpcg and eigh (using B as the metric) was within ~ 10^-8 of 1 (not surprising as this is the square of the first number). I ran a quick check where I dropped the imaginary parts of A and B and re-ran this analysis and say the errors fall to respectively ~ 10^-6 and 10^-12, so the algorithm gets less accurate results for complex numbers than for real numbers, though I don't know enough about how it works to speculate on which this would be. Anyway, I don't know much about lobpcg so there might be some issue that I've missed in simply adding ".conj()" everywhere a ".T" appeared to fix the dot products. I do think it would be very nice, though, to be able to use lobpcg for complex matrices, and so I would be willing to submit a patch towards this end. Thoughts? Cheers, Greg

3 2

Bundling Boost?
by Pauli Virtanen 12 Oct '12

12 Oct '12

Hi, I'd like to consider replacing some of the function implementations in scipy.special with versions from the C++ Boost library, cf. http://www.boost.org/doc/libs/1_51_0/libs/math/doc/sf_and_dist/html/index.h… The reason is that these implementations appear to be made with more care than what we use currently, and behave more correctly in corner cases (cf. e.g. ticket #1740). To minimize work, it would be useful just to use the Boost functions direcly, rather than doing manual C transcriptions. The drawback here is that the portion of Boost library required weights about 8 MB of source code, and we would most likely like to bundle it, as it is not really a standard part of many installations. This does not reflect much on the compiled binary size, however. I'm not 100 % certain about the compiler support. Perhaps C++ is already mature enough to work across the platforms we care about. I'm not aware of many good BSD-compatible floating-point special function libraries, so if you know others, or would be opposed to bundling Boost, please chime up! -- Pauli Virtanen

13 31

scipy.org will be down today for about 30 minutes
by Ognen Duzlevski 12 Oct '12

12 Oct '12

Hello, The scipy.org machine will be down today for about 30 minutes or so. I will be adding an extra drive to the machine since it has been slowly running out of space. Thanks! Ognen

1 0

Scipy docstrings
by Cera, Tim 10 Oct '12

10 Oct '12

A whole bunch of conflicts and merges showed up on the docstring editor at http://docs.scipy.org/scipy/merge/ Can't 'Accept Merges' because all of the merges that I looked at did not follow the docstring standard. There isn't a way to 'Reject Merges', so what I have done in the past is accept, then revert. Could this be handled in some other way? Kindest regards, Tim

4 5