Mailman 3 February 2014 - SciPy-Dev

Re: [SciPy-Dev] removing umfpack wrapper
by Robert Cimrman 24 Nov '14

24 Nov '14

On 11/12/2013 11:22 PM, Ralf Gommers wrote: > On Tue, Nov 12, 2013 at 10:20 AM, Robert Cimrman <cimrman3(a)ntc.zcu.cz>wrote: > >> Hi Ralf, >> >> On 11/10/2013 04:36 PM, Ralf Gommers wrote: >>> Hi all, >>> >>> This is long overdue, but I think it's time to remove the UMFPACK wrapper >>> before the next release. I'd like some feedback on how to do that and to >>> what package (if any) to point existing users. >>> >>> As for why to remove the wrapper, see: >>> https://github.com/scipy/scipy/issues/3002 >>> http://permalink.gmane.org/gmane.comp.python.scientific.user/20451 >>> Short summary: UMFPACK changed from a BSD-compatible license to GPL at >> some >>> point. >>> >>> The deprecation warning in sparse.linalg has been referring people to >>> scikits.umfpack until now, however that package has disappeared >> completely >>> as far as I can tell. I suspect it was in the old scikits svn repo and >> was >>> never moved before that was killed. The alternatives seems to be >> Pysparse ( >> >> I missed that the scikit was lost in transition... >> >> How much work would be to move it to the new site - move from SVN to git, >> and? >> I would be willing to do it, if there is interest. >> > > Hi Robert, I think the moving wouldn't be a lot of work (assuming svn > access can still be arranged). If you'd revive the scikit it's not a > one-time effort though - it's only useful if the code is being maintained > imho. That can still be low effort perhaps, but with reviewing some PRs, > maintenance releases, putting it on pypi, etc. I'd expect at least an hour > a week or so. > > Are there scikits.umfpack users now? Would it make sense to salvage > whatever is useful from the scikit and contribute it to pysparse instead? We are using it in sfepy as the default direct solver. It seems to me that the pysparse interface requires the input matrix to be in the LL (linked list) format, which is unfortunate for us, as we use CSR. The scipy (and former scikit) umfpack wrappers have used CSR, so no copies were necessary. (After all, I have created the original scipy wrappers to be used from sfepy in the first place...) So in case others are not interested in having the scikit, I can see two possible solutions from my perspective: either enhance pysparse interface to allow CSR as well, or move the wrappers to sfepy (which I maintain anyway, but swig -> cython conversion would be needed). Not sure yet which solution I prefer. Cheers, r. > Cheers, > Ralf >> >> Cheers, >> r. >> >>> http://pysparse.sourceforge.net) and PyUBlasExt ( >>> https://pypi.python.org/pypi/PyUblasExt). Scikits.sparse is dead(-ish) >> and >>> doesn't wrap umfpack it looks like. PyUBlasExt has something that looks >>> like a wrapper but for C++ (?). So can we refer people to pysparse? >>> >>> Other question: do we call Pysparse instead of scikits.umfpack when >>> available? And if so, by default or not? Or do we rip out everything that >>> sits under useUmfpack=True? >>> >>> Cheers, >>> Ralf

2 5

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 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

Levenberg-Marquardt Implementation
by Nico Del Piano 19 Jun '14

19 Jun '14

Hi all, I am Nicolas Del Piano, and currently studying the last year of computer science career. I was searching for a GSoC project that involves mathematical concepts, and I found one very interesting. I have experience on python programming and I studied calculus, so I have a background of the problem context. I am interested on the implementation and testing of Levenberg-Marquardt / trust region nonlinear least squares algorithm. Here are some useful links, that I have searched to having a reference about the problem: http://ananth.in/docs/lmtut.pdf (Introduction) www.cs.nyu.edu/~roweis/notes/lm.pdf (Optimization) http://scribblethink.org/Computer/Javanumeric/index.html (java implementation) I would be glad if there is an interested mentor to discuss the issues, talk about the problem and its possible implementations, and provide me some support and guide. Thanks! Regards.

10 21

conv2 and normxcorr2
by awebster＠falsecolour.com 12 May '14

12 May '14

I noticed a couple of popular matlab functions - conv2 and normxcorr2 were not present in the scipy.signal packages. I would like to submit them for addition. Can anyone point me on instructions on how to write such a thing? Below are examples. # 2D convolution using the convolution's FFT property def conv2(a,b): ma,na = a.shape mb,nb = b.shape return fft.ifft2(fft.fft2(a,[2*ma-1,2*na-1])*fft.fft2(b,[2*mb-1,2*nb-1])) # compute a normalized 2D cross correlation using convolutions # this will give the same output as matlab, albeit in row-major order def normxcorr2(b,a): c = conv2(a,flipud(fliplr(b))) a = conv2(a**2, ones(b.shape)) b = sum(b.flatten()**2) c = c/sqrt(a*b) return c -- Aaron Webster

5 8

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

Hankel transforms, again
by Tom Grydeland 23 Mar '14

23 Mar '14

Hi developers, This is a repost of a message from December 2008 which gave no useful answers. Since then, I’ve had 4-5 requests for the code from people who had a need for it. It’s not a massive demand, but enough that perhaps you’ll consider my offer again. Since the previous posting, I’ve also included alternative filters thanks to Fan-Nian Kong that are shorter and more accurate when the function makes significant changes in more limited intervals. I’m not including the code (since it is mostly thousands of lines of tables), but I will provide the files to anyone who’s interested. Cheers, Tom —— original message below ——— When I recently needed a Hankel transform I was unable to find one in Scipy. I found one in MATLAB though[1], written by Prof. Brian Borchers at New Mexico Tech. The code is based on a previous FORTRAN implementation by W. L. Anderson [2], and the MATLAB code is not marked with any copyright statements. Hankel transforms of the first and second order can be computed through digital filtering. I have rewritten the code in Python/numpy, complete with tests and acknowledgements of the origins, and my employer has agreed that I can donate the code to Scipy. I believe this should be of interest. Hankel transforms arise often in acoustics and other systems with cylinder symmetry. If you want it, I would like a suggestion where it belongs (with other integral transforms, probably) and how I should shape the tests to make them suitable for inclusion. The tests I currently have use Hankel transforms of five different functions with known analytical transforms and compares the transformed values to the numerically evaluated analytical expressions. Currently plots are generated, but for automated testing I suppose something else would be better. Pointing me at an example is sufficient. [1] http://infohost.nmt.edu/~borchers/hankel.html [2] Anderson, W. L., 1979, Computer Program Numerical Integration of Related Hankel Transforms of Orders 0 and 1 by Adaptive Digital Filtering. Geophysics, 44(7):1287-1305. Best regards, -- Tom Grydeland <Tom.Grydeland(a)(norut.no|gmail.com)>

4 13

API suggestions wanted for an enhancement to scipy.signal.filtfilt
by Warren Weckesser 10 Mar '14

10 Mar '14

Hey all, I'm adding an option to `scipy.signal.filtfilt` that uses Gustaffson's method [1] to handle the edges of the data. In this method, initial conditions for the forward and backward passes of `lfilter` are chosen such that applying the filter first in the forward direction and then backward gives the same result as applying the filter backward and then forward. There is no padding applied to the edges. Gustaffson's method has one optional parameter. It is the estimate of the length of the impulse response of the filter (i.e. anything after this length is supposed to be negligible and is ignored). If it is not given, no truncation of the impulse response is done. The current signature of `filtfilt` is def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None) (See http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.h… ) The arguments `padtype` and `padlen` control the type ('odd', 'even', 'constant' or None) and length of the padding. Any suggestions for modifying the signature in a backwards-compatible way? Here are a few options I've considered: (1) Specify the use of Gustaffson's method with `padtype='gust'`, and specify the estimate of the impulse response length using `padlen`. (I don't like this version--there is no padding performed by Gustaffson's method; using `padlen` for the impulse response length is just wrong.) (2) Specify the use of Gustaffson's method with `padtype='gust'`, and specify the estimate of the impulse response with a new argument `irlen`. (A bit better than (1); I could live with using `padtype` to specify the method, even though it isn't actually padding the data.) New signature: def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, irlen=None) (3) A new argument `method` specifies the method. It accepts either `'gust'` or `'pad'`. If the method is `'gust'`, the argument `irlen` specifies the impulse response length (and `padtype` and `padlen` are ignored). If the method is `'pad'`, `padtype` specifies the type of padding, `padlen` specifies the padding length (and `irlen` is ignored). The new signature would be: def filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad', irlen=None) (4) Don't touch `filtfilt`. Create a new function with the signature: def filtfilt_gust(b, a, x, axis=-1, irlen=None) Any suggestions? Any other APIs that would be preferable? Warren [1] F. Gustaffson. Determining the initial states in forward-backward filtering. Transactions on Signal Processing, 46(4):988-992, 1996.

2 2

0.13.3 and 0.14.0 releases
by Ralf Gommers 01 Mar '14

01 Mar '14

Hi all, It looks to me like we should do an 0.13.3 bugfix release soon to fix these two issues: - another memory leak in ndimage.label: https://github.com/scipy/scipy/issues/3148 - we broke weave.inline with Visual Studio: https://github.com/scipy/scipy/issues/3216 I propose to make this release within a week. For the 0.14.0 release I propose to branch this around Feb 23rd. That gives us a month to work through a decent part of the backlog of PRs. (plus I'm on holiday until the 15th, so earlier wouldn't work for me). Does that schedule work for everyone? Any more fixes that have to go into 0.13.3? Cheers, Ralf

5 8

GSoC: Cython-rewrite and improvement for scipy.cluster
by Richard Tsai 27 Feb '14

27 Feb '14

Hi all, I am Richard Tsai, and currently a second year undergraduate student in Computer Science at Sun Yat-sen University. I wish to take part in this year's GSoC. I'm learning machine learning with scipy/sklearn now and I've ever contributed some code to SciPy since last term. I've read Ralf's [Roadmap to Scipy 1.0][1] and I'm interested in the `cluster` part. I want to help finish the cython-rewrite work and make some improvement for it as my GSoC project. I noticed that there's a `_vq_rewrite.pyx` in `scipy/cluster` but I think it still needs furthor work. I want to start with some issues related to `cluster` as a warm-up, and then try to re-implement the `cluster.vq` module in cython first and try to do some optimizations. I'm familiar with it since I've ever done a little SNS text mining research with it with my classmates in a contest. As for the `cluster.hierarchy` module, I do not know a lot about hierarchical clustering for I haven't used it in practice. I may start with reading some papers and writing some examples for the documents. Then I will start the cython-rewrite for the `hierarchy` module. Finally, I plan to make some enhancements for the package. Maybe to automatically determine the number of clusters with Elbow Method? I haven't had a detailed plan yet. Since this idea is not listed on the ideas page, I don't know if it is suitable to be a GSoC project. If you have any suggestions, please let me know. I'd appreciate it if you can provide any guidance/opinions/suggestions. Regards, Richard [1]: https://github.com/rgommers/scipy/blob/roadmap/doc/ROADMAP.rst.txt

2 1