Mailman 3 September 2013 - SciPy-Dev

Linalg: function for diagonalize hermitian commuting matrices simultaneously
by Lars-Hendrik Frahm 07 Oct '14

07 Oct '14

Dear SciPy Developers, I like to contribute a function for simultaneously diagonalize a set of hermitian commuting matrices. Especially for matrices with degenerate eigenvalues. Is there need for something like that in SciPy ? I think it could be helpful for someone else and I didn't find it in recent SciPy versions or on the mailing list. Hope to hear from you.

3 3

Bootstrap confidence limits code
by Constantine Evans 04 Oct '14

04 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

1.0 roadmap: weave
by Arnd Baecker 06 May '14

06 May '14

Hi, > On Sat, Sep 21, 2013 at 8:54 PM, Ralf Gommers <ralf.gommers(a)gmail.com> wrote: [...] > weave > ````` > This is the only module that was not ported to Python 3. Effectively it's > deprecated (not recommended to use for new code). In the future it should be > removed from scipy (can be made into a separate module). Let me just mention that I am worried about the removal of weave: a substantial amount of code developed in our group over the years uses weave. Porting this to cython seems not trivial, in particular, because many of the original authors (diploma students, PhD students, ...) have left academia by now. Well, as porting weave to Python 3 is way out of my capabilities, this is just a comment in the hope that someone could give this a try (if at some point testing would be needed, I could try to help, of course). Best, Arnd

9 12

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

Vector Strength function
by Todd 21 Nov '13

21 Nov '13

I am interested in implementing a function for scipy. The function is called "vector strength". It is basically a measure of how reliably a set of events occur at a particular phase. It was originally developed for neuroscience research, to determine how well a set of neural events sync up with a periodic stimulus like a sound waveform. However, it is useful for determining how periodic a supposedly periodic set of events really are, for example: 1. Determining whether crime is really more common during a full moon and by how much 2. Determining how concentrated visitors to a coffee shop are during rush hour 3. Determining exactly how concentrated hurricanes are during hurricane season My thinking is that this could be implemented in stages: First would be a Numpy function that would add a set of vectors in polar coordinates. Given a number of magnitude/angle pairs it would provide a summed magnitude/angle pair. This would probably be combined with a cartesian<->polar conversion functions. Making use of this function would be a scipy function that would actually implement the vector strength calculation. This is done by treating each event as a unit vector with a phase, then taking the average of the vectors. If all events have the same phase, the result will have an amplitude of 1. If they all have a different phases, the result will have an amplitude of 0. It may even be worth having a dedicated polar dtype, although that may be too much. What does everyone think of this proposal?

4 12

Quadpack conversion to Fortran
by Brian Lee Newsom 09 Oct '13

09 Oct '13

Hello, My name is Brian Newsom and I am just beginning to use/work on the scipy library, specifically the quadpack within the integration pack. Because of the recursion and wrapping methods used now, running nquad (scipy 0.13) function is slow-the fortran and python are nested in a way that requires switching back and forth between the two often. Ideally this inefficiency could be overcome by translating the main interface (quadpack.py) into fortran and just having a minimal python wrapper that is not used in the recursion. First off, has anyone attempted (or is currently attempting) this? Are there any major difficulties I should be wary of or is this unrealistic for any reason? Additionally, I was curious if there are any computational advantages to using C or Fortran for this purpose within this library, either is a possibility in the case that one is strongly preferred. Otherwise, I am sure I will be full of questions in the near future with specific regards to the building and setup of these files within the scipy library. Thank you, Brian Newsom Brian.Newsom(a)colorado.edu

3 5

1.8.0rc1
by Charles R Harris 02 Oct '13

02 Oct '13

Hi All, NumPy 1.8.0rc1 is up now on sourceforge<http://sourceforge.net/projects/numpy/files/NumPy/1.8.0rc1/>.The binary builds are included except for Python 3.3 on windows, which will arrive later. Many thanks to Ralf for the binaries, and to those who found and fixed the bugs in the last beta. Any remaining bugs are all my fault ;) I hope this will be the last release before final, so please test it thoroughly. Chuck

6 9

ImportError: cannot import name oldnumeric
by xavier.gnata＠gmail.com 29 Sep '13

29 Sep '13

Hi, I have compiled numpy1.9.0.dev-e051ff following http://osdf.github.io/blog/numpyscipy-with-openblas-for-ubuntu-1204-second-… to unable Openblas. It looks ok (run=5277 errors=0 failures=0 and Openblas is used) but now scipy git complains like this: "python setup.py build Cythonizing sources scipy/interpolate/interpnd.pyx has not changed scipy/ndimage/src/_ni_label.pyx has not changed scipy/stats/vonmises_cython.pyx has not changed scipy/stats/_rank.pyx has not changed scipy/signal/_spectral.pyx has not changed scipy/special/_ufuncs_cxx.pyx has not changed scipy/special/_ufuncs.pyx has not changed scipy/cluster/_vq_rewrite.pyx has not changed scipy/sparse/csgraph/_shortest_path.pyx has not changed scipy/sparse/csgraph/_traversal.pyx has not changed scipy/sparse/csgraph/_min_spanning_tree.pyx has not changed scipy/sparse/csgraph/_tools.pyx has not changed scipy/io/matlab/mio5_utils.pyx has not changed scipy/io/matlab/mio_utils.pyx has not changed scipy/io/matlab/streams.pyx has not changed scipy/spatial/ckdtree.pyx has not changed scipy/spatial/qhull.pyx has not changed /usr/lib/python2.7/distutils/dist.py:267: UserWarning: Unknown distribution option: 'test_suite' warnings.warn(msg) blas_opt_info: blas_mkl_info: Disabled blas_mkl_info: (MKL is None) Disabled blas_mkl_info: (MKL is None) libraries mkl,vml,guide not found in [] NOT AVAILABLE openblas_info: FOUND: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] language = f77 FOUND: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] language = f77 Traceback (most recent call last): File "setup.py", line 230, in <module> setup_package() File "setup.py", line 227, in setup_package setup(**metadata) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/core.py", line 135, in setup config = configuration() File "setup.py", line 170, in configuration config.add_subpackage('scipy') File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 966, in add_subpackage caller_level = 2) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 935, in get_subpackage caller_level = caller_level + 1) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 872, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/setup.py", line 15, in configuration config.add_subpackage('lib') File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 966, in add_subpackage caller_level = 2) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 935, in get_subpackage caller_level = caller_level + 1) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 872, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/lib/setup.py", line 9, in configuration config.add_subpackage('blas') File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 966, in add_subpackage caller_level = 2) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 935, in get_subpackage caller_level = caller_level + 1) File "/usr/local/lib/python2.7/dist-packages/numpy/distutils/misc_util.py", line 872, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/lib/blas/setup.py", line 22, in configuration from scipy._build_utils import get_g77_abi_wrappers File "/usr/local/src/scipy/scipy/__init__.py", line 77, in <module> from numpy import oldnumeric ImportError: cannot import name oldnumeric " What's going on with oldnumeric? xavier

2 2

Scipy 1.0 roadmap - wavelets
by Daπid 28 Sep '13

28 Sep '13

On 21 September 2013 20:57, Ralf Gommers <ralf.gommers(a)gmail.com> wrote: > *Wavelets*: add proper wavelets, including discrete wavelet transform. > What's > there now doesn't make much sense. > There is a package on Wavelets under MIT license, but it looks abandoned for more than a year. Perhaps it could be reintegrated into Scipy, perhaps it would be better to create a scikit-wavelets while we get something of enough quality. http://www.pybytes.com/pywavelets I can step in, if someone helps and guides me. /David.

3 4

Re: [SciPy-Dev] [Numpy-discussion] Baffling error: ndarray += csc_matrix -> "ValueError: setting an array element with a sequence"
by Nathaniel Smith 27 Sep '13

27 Sep '13

On Fri, Sep 27, 2013 at 8:27 PM, Pauli Virtanen <pav(a)iki.fi> wrote: > 27.09.2013 22:15, Nathaniel Smith kirjoitti: > [clip] >> 3) The issue of how to make an in-place like ndarray += sparse >> continue to work in the brave new __numpy_ufunc__ world. >> >> For this last issue, I think we disagree. It seems to me that the >> right answer is that csc_matrix.__numpy_ufunc__ needs to step up and >> start supporting out=! If I have a large dense ndarray and I try to += >> a sparse array to it, this operation should take no temporary memory >> and nnz time. Right now it sounds like it actually copies the large >> dense ndarray, which takes time and space proportional to its size. >> AFAICT the only way to avoid that is for scipy.sparse to implement >> out=. It shouldn't be that hard...? > > Sure, scipy.sparse can easily support also the output argument. Great! I guess solving this somehow will be release-critical, to avoid a regression in this case when __numpy_ufunc__ gets released. If the solution should be in scipy, I guess we should file the bug there? > But I still maintain that the implementation of __iadd__ in Numpy is > wrong. Oh yeah totally. > What it does now is: > > def __iadd__(self, other): > return np.add(self, other, out=self) > > But since we know this raises a TypeError if the input is of a type that > cannot be dealt with, it should be: > > def __iadd__(self, other): > try: > return np.add(self, other, out=self) > except TypeError: > return NotImplemented > > Of course, it's written in C so it's a bit more painful to write this. > > I think this will have little performance impact, since the check would > be only a NULL check in the inplace op methods + subsequent handling. I > can take a look at some point at this... I'm a little uncertain about the "swallow all TypeErrors" aspect of this -- e.g. this could have really weird effects for object arrays, where ufuncs may raise arbitrary user exceptions. One possibility in the long run is to just say, if you want to override ndarray __iadd__ or whatever, then you have to use __numpy_ufunc__. Not really much point in having *two* sets of implementations of the NotImplemented dance for the same operation. -n

2 1