Mailman 3 July 2011 - SciPy-User

[ANN] IPython 0.11 is officially out
by Fernando Perez 03 Aug '11

03 Aug '11

Hi all, on behalf of the IPython development team, I'm thrilled to announce, after more than two years of development work, the official release of IPython 0.11. This release brings a long list of improvements and new features (along with hopefully few new bugs). We have completely refactored IPython, making it a much more friendly project to participate in by having better separated and organized internals. We hope you will not only use the new tools and libraries, but also join us with new ideas and development. After this very long development effort, we hope to make a few stabilization releases at a quicker pace, where we iron out the kinks in the new APIs and complete some remaining internal cleanup work. We will then make a (long awaited) IPython 1.0 release with these stable APIs. *Downloads* Download links and instructions are at: http://ipython.org/download.html And IPython is also on PyPI: http://pypi.python.org/pypi/ipython Those contain a built version of the HTML docs; if you want pure source downloads with no docs, those are available on github: Tarball: https://github.com/ipython/ipython/tarball/rel-0.11 Zipball: https://github.com/ipython/ipython/zipball/rel-0.11 * Features * Here is a quick listing of the major new features: - Standalone Qt console - High-level parallel computing with ZeroMQ - New model for GUI/plotting support in the terminal - A two-process architecture - Fully refactored internal project structure - Vim integration - Integration into Microsoft Visual Studio - Improved unicode support - Python 3 support - New profile model - SQLite storage for history - New configuration system - Pasting of code with prompts And many more... We closed over 500 tickets, merged over 200 pull requests, and more than 60 people contributed over 2200 commits for the final release. Please see our release notes for the full details on everything about this release: https://github.com/ipython/ipython/zipball/rel-0.11 As usual, if you find any problem, please file a ticket --or even better, a pull request fixing it-- on our github issues site (https://github.com/ipython/ipython/issues/). Many thanks to all who contributed! Fernando, on behalf of the IPython development team. http://ipython.org

1 1

deconvolution of 1-D signals
by Ralf Gommers 01 Aug '11

01 Aug '11

Hi, For a measured signal that is the convolution of a real signal with a response function, plus measurement noise on top, I want to recover the real signal. Since I know what the response function is and the noise is high-frequency compared to the real signal, a straightforward approach is to smooth the measured signal (or fit a spline to it), then remove the response function by deconvolution. See example code below. Can anyone point me towards code that does the deconvolution efficiently? Perhaps signal.deconvolve would do the trick, but I can't seem to make it work (except for directly on the output of np.convolve(y, window, mode='valid')). Thanks, Ralf import numpy as np from scipy import interpolate, signal import matplotlib.pyplot as plt # Real signal x = np.linspace(0, 10, num=201) y = np.sin(x + np.pi/5) # Noisy signal mode = 'valid' window_len = 11. window = np.ones(window_len) / window_len y_meas = np.convolve(y, window, mode=mode) y_meas += 0.2 * np.random.rand(y_meas.size) - 0.1 if mode == 'full': xstep = x[1] - x[0] x_meas = np.concatenate([ \ np.linspace(x[0] - window_len//2 * xstep, x[0] - xstep, num=window_len//2), x, np.linspace(x[-1] + xstep, x[-1] + window_len//2 * xstep, num=window_len//2)]) elif mode == 'valid': x_meas = x[window_len//2:-window_len//2+1] elif mode == 'same': x_meas = x # Approximating spline xs = np.linspace(0, 10, num=500) knots = np.array([1, 3, 5, 7, 9]) tck = interpolate.splrep(x_meas, y_meas, s=0, k=3, t=knots, task=-1) ys = interpolate.splev(xs, tck, der=0) # Find (low-frequency part of) original signal by deconvolution of smoothed # measured signal and known window. y_deconv = signal.deconvolve(ys, window)[0] #FIXME # Plot all signals fig = plt.figure() ax = fig.add_subplot(111) ax.plot(x, y, 'b-', label="Original signal") ax.plot(x_meas, y_meas, 'r-', label="Measured, noisy signal") ax.plot(xs, ys, 'g-', label="Approximating spline") ax.plot(xs[window.size//2-1:-window.size//2], y_deconv, 'k-', label="signal.deconvolve result") ax.set_ylim([-1.2, 2]) ax.legend() plt.show()

9 18

Re: [SciPy-User] integral of oscillatory functions
by Juan Fiol 31 Jul '11

31 Jul '11

Thanks for all the answers. They were very helpful. In general, seems that there is no other way that to tackle the specific problem with some analytical work. I'l look more deeply into that, and let you know if something interesting pops out. Joshua, thanks for the link. I think that the method at what the document refers is the one I've mentioned as "asymptotic". I'd already skimmed over most of the papers cited there but as I said only the surface. I think that this document has an amenable enough presentation that will make it useful to give a try. Adam: I'll look further into the asymptotic methods. The problem is that the oscillatory part is complicated enough to make it painful to go through the analytical work. I do not plan to do the work in scipy/numpy. I am trying to solve the problem in python but need then to adapt the solution to one of my fortran programs. Charles: I apologize for not being specific enough. The integrand itself is quite complicated. Moreover its form may change. I am attaching a short pdf with my current definitions, but in python would be something as # Auxiliary functions and definitions # These values will be changing w0= 0.4 dw= 0.1 w=np.array([w0,w0+dw, w0-dw]) A=3. B=4. Omega= B - (A*(a/np.square(w))).sum() def h(t): return ((a/w)*(np.cos(w*t))).sum() + C def G(t): return (k_A*(a/np.square(w))*np.sin(w*t)).sum() # function to integrate def integrand(t): return np.exp(1j* G(t)) * np.exp(1j* Omega*t) * f(t) --- On Sat, 7/30/11, Charles R Harris <charlesr.harris(a)gmail.com> wrote: From: Charles R Harris <charlesr.harris(a)gmail.com> Subject: Re: [SciPy-User] integral of oscillatory functions To: fiolj(a)yahoo.com, "SciPy Users List" <scipy-user(a)scipy.org> Date: Saturday, July 30, 2011, 1:10 AM On Fri, Jul 29, 2011 at 4:36 PM, Juan Fiol <fiolj(a)yahoo.com> wrote: Hi, I have to integrate a *highly* oscilatory function. I've been looking in the literature and found that there some "asymptotic methods" (that work better when oscillations are stronger and cancel most of the integrand), some methods derived from filon's method and other called Levin method. I've had taken quick looks into several mathematical papers on the subject but It will probably take me more than one month (and may be much more) to understand the subject and put it into a routine. Does anybody know if there is anything of the sort implemented in scipy? Otherwise, I would appreciate if I get advice for a more "practical" place where to look. The integrand is not strictly of the form f(t) e^(iwt). Any help would be welcome. Thanks How oscillatory is *highly* oscillatory? Does the function have any particular form? Where does it come from? Chuck

3 2

recommendation for saving data
by Brian Blais 31 Jul '11

31 Jul '11

Hello, I was wondering if there are any recommendations for formats for saving scientific data. I am running a simulation, which has many somewhat-indepedent parts which have their own internal state and parameters. I've been using pickle (gzipped) to save the entire object (which contains subobjects, etc...), but it is getting too unwieldy and I think it is time to look for a more robust solution. Ideally I'd like to have something where I can call a save method on the simulation object, and it will call the save methods on all the children, on down the line all saving into one file. It'd also be nice if it were cross-platform, and I could depend on the files being readable into the future for a while. Are there any good standards for this? What do you use for saving scientific data? thank you, Brian Blais -- Brian Blais bblais(a)bryant.edu http://web.bryant.edu/~bblais http://bblais.blogspot.com/

2 1

creating a view of an array
by mathieu lacage 30 Jul '11

30 Jul '11

hi, Let's say I have a big large array: a = numpy.empty((10000*1000,10)) and I want to create a view of that array to be able to process a subset of its data without making a copy. One column: b = a[:,1] Two adjacent columns: b = a[:,1:2] Can I do the same (no memory allocation) for two columns that are not adjacent ? i.e. if I want to create an array that is a view for the 1st and 3rd columns. Mathieu

4 4

generic_flapack.pyf and geqp3
by Collin Stocks 30 Jul '11

30 Jul '11

Hi, all. I am planning to try to add functionality to scipy.linalg.qr(), specifically to allow qr decomposition with pivoting. However, I have almost no knowledge of how to wrap the function in scipy/linalg/ generic_flapack.pyf. Could somebody please point me in the correct direction? Thanks, Collin Stocks

2 6

Sparse eigensystem instability on OS X Lion
by Chris Laumann 30 Jul '11

30 Jul '11

Hi all- I've recently upgraded to Lion and, because it broke my python setup, reinstalled scipy from two sources -- Enthought's full 64 bit latest distro (Scipy 0.9.0 and Numpy 1.6.0) and Chris Fonnesbeck's Scipy Superpack (with the latest dev versions from a few weeks ago). I'm now finding a bunch of numerical instability in symmetric sparse eigensystem solutions in scipy.sparse.linalg.eigsh. These didn't exist before the upgrade when I was using Snow Leopard and scipy 0.8.0 (eigen_symmetric instead of eigsh). The smallest example I've isolated is: data = array([ 16., 18., -4., -2., -4., 25., -2., 4., -2., 17., 16., 18.]) indices = array([0, 1, 2, 6, 1, 2, 4, 3, 2, 4, 5, 6], dtype='int32') indptr = array([ 0, 1, 4, 7, 8, 10, 11, 12], dtype='int32') ham1 = scipy.sparse.csr_matrix((data, indices, indptr)) e1 = np.linalg.eigvalsh(ham1.todense()) e2 = scipy.sparse.linalg.eigsh(ham1, k=6, which='SA', return_eigenvectors = False) e2.sort() This leads to (using the superpack -- the results from Enthought are similarly noisy but different): >>> e1 array([ 4. , 15.57851688, 16. , 16. , 17.27739936, 18. , 27.14408376]) >>> e2 array([ 4.03171393, 15.47517729, 15.98093481, 16. , 16.76781718, 18.79560773]) These are big errors for a 7x7 matrix, especially for the minimum eigenvalue 4, which is well separated from the rest of the spectrum. Can anybody help? I believe the Superpack _arpack.so is linking against the apple Accelerate framework for LAPACK while I'm not sure about Enthoughts _arpack.so. As far as I can tell from a little otool snooping, Enthought seems to roll its own LAPACK functionality. In any event, neither seems to work. Thanks, Chris

2 1

integral of oscillatory functions
by Juan Fiol 29 Jul '11

29 Jul '11

Hi, I have to integrate a *highly* oscilatory function. I've been looking in the literature and found that there some "asymptotic methods" (that work better when oscillations are stronger and cancel most of the integrand), some methods derived from filon's method and other called Levin method. I've had taken quick looks into several mathematical papers on the subject but It will probably take me more than one month (and may be much more) to understand the subject and put it into a routine. Does anybody know if there is anything of the sort implemented in scipy? Otherwise, I would appreciate if I get advice for a more "practical" place where to look. The integrand is not strictly of the form f(t) e^(iwt). Any help would be welcome. Thanks Juan

4 3

DISREGARD: Sparse eigensystem instability on OS X Lion
by Chris Laumann 29 Jul '11

29 Jul '11

Hey everybody- The matrix I sent out wasn't hermitian, hence the instability. Disregard the previous post.. Oops, Chris ---------------- Hi all- I've recently upgraded to Lion and, because it broke my python setup, reinstalled scipy from two sources -- Enthought's full 64 bit latest distro (Scipy 0.9.0 and Numpy 1.6.0) and Chris Fonnesbeck's Scipy Superpack (with the latest dev versions from a few weeks ago). I'm now finding a bunch of numerical instability in symmetric sparse eigensystem solutions in scipy.sparse.linalg.eigsh. These didn't exist before the upgrade when I was using Snow Leopard and scipy 0.8.0 (eigen_symmetric instead of eigsh). The smallest example I've isolated is: data = array([ 16., 18., -4., -2., -4., 25., -2., 4., -2., 17., 16., 18.]) indices = array([0, 1, 2, 6, 1, 2, 4, 3, 2, 4, 5, 6], dtype='int32') indptr = array([ 0, 1, 4, 7, 8, 10, 11, 12], dtype='int32') ham1 = scipy.sparse.csr_matrix((data, indices, indptr)) e1 = np.linalg.eigvalsh(ham1.todense()) e2 = scipy.sparse.linalg.eigsh(ham1, k=6, which='SA', return_eigenvectors = False) e2.sort() This leads to (using the superpack -- the results from Enthought are similarly noisy but different): >>> e1 array([ 4. , 15.57851688, 16. , 16. , 17.27739936, 18. , 27.14408376]) >>> e2 array([ 4.03171393, 15.47517729, 15.98093481, 16. , 16.76781718, 18.79560773]) These are big errors for a 7x7 matrix, especially for the minimum eigenvalue 4, which is well separated from the rest of the spectrum. Can anybody help? I believe the Superpack _arpack.so is linking against the apple Accelerate framework for LAPACK while I'm not sure about Enthoughts _arpack.so. As far as I can tell from a little otool snooping, Enthought seems to roll its own LAPACK functionality. In any event, neither seems to work. Thanks, Chris

1 0

What happened with chisquare_contingency?
by Yang Zhang 29 Jul '11

29 Jul '11

Looked like a useful addition: http://mail.scipy.org/pipermail/scipy-dev/2010-June/014538.html I originally was rolling my own and found the chisquare ddof parameter to be confusing - Googling turned this up. Just wondering if this is still on its way into the codebase. Thanks!

4 4