Mailman 3 December 2008 - SciPy-Dev

fast gaussian filtering
by Sturla Molden Dec. 19, 2008

Dec. 19, 2008

Gaussian filtering is a very useful method for smoothing signals and images, as well as for kernel density estimation. I have previously used it to estimate the firing rates of neurons. I've noticed that SciPy implements Gaussian filtering (in the ndimage package) using convolution. It is implemented by loops in Python, not even using a vectorized numpy.convolve or an FFT. For the sake of speed, SciPy is using the worst thinkable implementation. To make it fast: Gaussian filtering can be implemented using recursion (Young & Vliet, 1995; Signal Processing, 44: 139-151). The recursive approximation is very accurate and does not introduce ringing. It is anti-causal (forward-backward) and has zero phase response. The filtering can be done in C using scipy's signal.lfilter function. I found a solution in plain C on the Matlab central. But there is really no benefit from doing this in C as lfilter, flipud and fliplr does all the hard work. An efficient Matlab version could have been constructed similarly, without resorting to C. http://www.mathworks.com/matlabcentral/fileexchange/5087 http://www.ph.tn.tudelft.nl/~lucas/publications/papersLJvV.html Note that Young & Vliet (1995) also provides recursive solutions for the first, second and third derivatives. They are equally easy to implement using lfilter. Regards, Sturla Molden from numpy import array, zeros, ones, flipud, fliplr from scipy.signal import lfilter from math import sqrt def __gausscoeff(s): if s < .5: raise ValueError, \ 'Sigma for Gaussian filter must be >0.5 samples' q = 0.98711*s - 0.96330 if s > 0.5 else 3.97156 \ - 4.14554*sqrt(1.0 - 0.26891*s) b = zeros(4) b[0] = 1.5785 + 2.44413*q + 1.4281*q**2 + 0.422205*q**3 b[1] = 2.44413*q + 2.85619*q**2 + 1.26661*q**3 b[2] = -(1.4281*q**2 + 1.26661*q**3) b[3] = 0.422205*q**3 B = 1.0 - ((b[1] + b[2] + b[3])/b[0]) # convert to a format compatible with lfilter's # difference equation B = array([B]) A = ones(4) A[1:] = -b[1:]/b[0] return B,A def Gaussian1D(signal, sigma, padding=0): n = signal.shape[0] tmp = zeros(n + padding) if tmp.shape[0] < 4: raise ValueError, \ 'Signal and padding too short' tmp[:n] = signal B,A = __gausscoeff(sigma) tmp = lfilter(B, A, tmp) tmp = tmp[::-1] tmp = lfilter(B, A, tmp) tmp = tmp[::-1] return tmp[:n] def Gaussian2D(image, sigma, padding=0): n,m = image.shape[0],image.shape[1] tmp = zeros((n + padding, m + padding)) if tmp.shape[0] < 4: raise ValueError, \ 'Image and padding too small' if tmp.shape[1] < 4: raise ValueError, \ 'Image and padding too small' B,A = __gausscoeff(sigma) tmp[:n,:m] = image tmp = lfilter(B, A, tmp, axis=0) tmp = flipud(tmp) tmp = lfilter(B, A, tmp, axis=0) tmp = flipud(tmp) tmp = lfilter(B, A, tmp, axis=1) tmp = fliplr(tmp) tmp = lfilter(B, A, tmp, axis=1) tmp = fliplr(tmp) return tmp[:n,:m]

2 3

new incremental statistics project
by Bradford Cross Dec. 19, 2008

Dec. 19, 2008

This is a new project I just released. I know it is C#, but some of the design and idioms would be nice in numpy/scipy for working with discrete event simulators, time series, and event stream processing. http://code.google.com/p/incremental-statistics/

1 0

hidden scipy.optimize.nnls
by josef.pktd＠gmail.com Dec. 18, 2008

Dec. 18, 2008

Is scipy.optimize.nnls hidden on purpose? It is not included in optimize.__init__.__all__ nor in info.py and is also not in the docs. Or should we add to optimize.__init__? from nnls import nnls Josef

3 3

Scipy Tutorial (and updating it)
by Pauli Virtanen Dec. 18, 2008

Dec. 18, 2008

Dear all, I've gone through a first pass of updating the old Scipy Tutorial to present day. The result is here (and sources in Scipy SVN trunk/doc/ source/tutorial): http://docs.scipy.org/doc/scipy/reference/tutorial/index.html It would be good if some brave souls went through it and helped to check if there are still some bugs hiding in the text. You can also fix any issues you find immediately here: http://docs.scipy.org/scipy/docs/scipy-docs/tutorial/index.rst/ provided you have an account on the site. Note that plots won't be visible on the wiki page -- you'll need to make sure that the stuff you have in the plot:: directives produces the appropriate plots manually. -- Pauli Virtanen

7 28

Re: [SciPy-dev] FFT docstrings (was: Scipy Tutorial (and updating it))
by jh＠physics.ucf.edu Dec. 18, 2008

Dec. 18, 2008

Tom Grydeland <Tom.Grydeland(a)(gmail.com)> writes: > I've started on the module page, along the direction you outlined > above, not complete and no examples quite yet, but if people think it > looks okay, I can try to trim down the specific routines a bit as > well. I like it! A lot! I added the helper routines and defined the time and frequency domains near the top (and took out similar text near the bottom), fixed some commas, etc. I think all it needs now is examples. One point: (I'm not sure it's worth putting these distinctions on the page, I'm just justifying a change I made.) Power is energy per unit time; the energy spectrum and power spectrum are colloquially the same. Both are proportional to the square of the amplitude spectrum. In areas where it matters, np.abs(A)**2 is the energy spectrum and varies with the signal length (longer signals have more energy). The power spectrum is the energy spectrum divided by the time span (in physical units) of the dataset, and is independent of the signal length. Both these must be adjusted by a constant to make the sampling rate not affect the output. The zero channel of np.fft(a) is the mean, not the energy. The zero channel of np.abs(A)**2 is the total energy or power, each to within different constants. --jh--

3 4

suggestion for Numpy documentation editor
by timmichelsen＠gmx-topmail.de Dec. 18, 2008

Dec. 18, 2008

Hello. I would like to ask whether the Numpy documentation editor could be implemented with a password reminder/resend option. Also, the follwoing links are dead: http://docs.scipy.org/numpy/Front%20Page/ => http://docs.scipy.org/doc/numpy.doc.reference.jargon http://docs.scipy.org/doc/numpy.doc.reference.glossary Thanks and regards, Timmie -- Psssst! Schon vom neuen GMX MultiMessenger gehört? Der kann`s mit allen: http://www.gmx.net/de/go/multimessenger

1 0

ANNOUNCE: EPD Py25 v4.1.30101_beta2 available for testing
by Dave Peterson Dec. 18, 2008

Dec. 18, 2008

Hello, The Enthought Python Distribution's (EPD) early access program website is now hosting the beta 2 build of the upcoming EPD Py25 v4.1.301 release. We would very much appreciate your assistance in making EPD as stable and reliable as possible! Please join us in our efforts by downloading an installer for Windows, Mac OS X, or RedHat EL versions 3, 4, and 5 from the following website: http://www.enthought.com/products/epdearlyaccess.php The release notes for the beta2 build are available here: https://svn.enthought.com/epd/wiki/Py25/4.1.301/Beta2 Please provide any comments, concerns, or bug reports via the EPD Trac instance at https://svn.enthought.com/epd or via e-mail to epd-support(a)enthought.com. -- Dave About EPD --------- The Enthought Python Distribution (EPD) is a "kitchen-sink-included" distribution of the Python™ Programming Language, including over 60 additional tools and libraries. The EPD bundle includes NumPy, SciPy, IPython, 2D and 3D visualization, database adapters, and a lot of other tools right out of the box. http://www.enthought.com/products/epd.php It is currently available as a single-click installer for Windows XP (x86), Mac OS X (a universal binary for OS X 10.4 and above), and RedHat 3, 4, and 5 (x86 and amd64). EPD is free for academic use. An annual subscription and installation support are available for individual commercial use. Various workgroup, departmental, and enterprise level subscription options with support and training are also available. Contact us for more information!

1 0

scipy.test() abruptly quits on scipy svn 0.7.0.dev5241
by Peter Norton Dec. 16, 2008

Dec. 16, 2008

I've had the same result for the last few checkouts. I've been building with numscons, and as far as I can tell no severe errors are coming up in the build. However I get this: ....zswap:n=4 ..zswap:n=3 .... **************************************************************** WARNING: clapack module is empty ----------- See scipy/INSTALL.txt for troubleshooting. Notes: * If atlas library is not found by numpy/distutils/system_info.py, then scipy uses flapack instead of clapack. **************************************************************** ...Result may be inaccurate, approximate err = 1.91669384639e-08 ........SSS.................................................................................................................................................................................................................................................................................................................................................................................................................................................ESSSSSSSSSSS. ** On entry to CLASCL parameter number 4 had an illegal value $ echo $? 0 This looks like it's bailing out, but it doesn't appear to be exiting with any errors. This was invoked as: $ python -c 'import scipy,sys; sys.exit(scipy.test())' The same error occurs whether I build with numscons or distutils (though numscons is a lot easier). Is this expected behavior? Thanks, -Peter

3 8

doc review & sphinx extension help
by timmichelsen＠gmx-topmail.de Dec. 16, 2008

Dec. 16, 2008

Dear Scipy developers, would you please enable my account "timmichelsen" for editing scipy and numpy documentation? I would like to add examples and review. Furthermore, I started with a bibtex extension for Sphinx. But I have problems with extracting citations from docstrings. Maybe the creators of the sphinx extension have a bright idea. I have a description of my approach here: * http://groups.google.com/group/sphinx-dev/browse_thread/thread/16ba98801e8e… I hope to paste an example within this week to http://www.bitbucket.org/birkenfeld/sphinx/issue/63/make-sphinx-read-bibtex… Thanks in advance, Timmie -- Psssst! Schon vom neuen GMX MultiMessenger gehört? Der kann`s mit allen: http://www.gmx.net/de/go/multimessenger

2 1

usage of "scipy", "Scipy", "SciPy"
by jh＠physics.ucf.edu Dec. 16, 2008

Dec. 16, 2008

Damian Eads wrote in a thread on scipy-user (ANN: SciKits Portal): > 2. Replace "Scipy" and "scipy" with "SciPy" when referring to the > SciPy tool suite. Use "scipy" (all lowercase) in a fixed-width font > when referring to the Python package name "scipy". David Cournapeau <david(a)ar.media.kyoto-u.ac.jp> responded: > I think it would be nice to have this noted somewhere for everything > around scipy, actually. +1, mostly. I don't think the fixed-width font is crucial; we don't consistently do this for "numpy" except in actual code examples and I for one find it jarring to read prose in manuals where you know perfectly well what is being discussed but fonts are changing all the time as though you're an idiot and can't remember. I do find it helpful to change fonts in syntax usages, etc., just not in paragraphs of prose when the context is clear (e.g., you're just talking about one item). Others will differ. So, I'd make the suggestion and leave it to the writer to decide on fonts appropriate to the document. Also, what about numpy/Numpy/NumPy? See below. We have a set of doc guidelines. I don't think we (yet) have a set of style guidelines for web sites, but we should. This could be its first point. ---------------------------------------------------------------------- Use "SciPy" when referring to the SciPy tool suite. When referring to the single Python package named "scipy", use "scipy" (all lowercase), preferably in a fixed-width font when using fonts to distinguish syntactic elements. At or very near first use, and in titles, etc., say "SciPy platform" when referring to the tool suite and "scipy package" when referring to the package alone, to ensure that readers not familiar with the distinction are not confused. When capitalization is required, "scipy" becomes "Scipy" and should be followed by "package" unless the context is already well established. Use "NumPy" when referring in general to the syntax, mechanics, and coding approach of the numpy package and/or the ndarray class, as used in other packages. When referring to the single Python package named "numpy", use "numpy" (all lowercase), preferably in a fixed-width font when using fonts to distinguish syntactic elements. Writers should clarify the context if there is ambiguity. When capitalization is required, "numpy" becomes "Numpy". Examples: The SciPy platform includes the numpy, scipy, and matplotlib packages. The SciPy universe is extended by the stand-alone packages called SciKits, and many other packages distributed separately from the SciKits. The fnord package's routines follow NumPy style. They use ndarrays and many accept "axis=" and other NumPy keywords. The fnord package is based on numpy. It also uses matplotlib, but it does not import the scipy package. (Title) The SciPy Platform User Manual (Title) The Scipy Package User Manual Incorrect: (Title) The SciPy User Manual (requires new reader to know our nuances) (Title) Implementation of the SciPy "axis=" Keyword (should be NumPy) In the NumPy package,... (should be numpy) In numpy/scipy, this IDL expression shortens to... (should be NumPy for generic expression syntax, SciPy if specific routines outside of numpy are involved, or numpy or scipy if routines from those packages are involved) ----------------------------------------------------------------------- Thoughts? I'm open to calling the syntactic elements "SciPy" rather than "NumPy", and never allowing "NumPy" or allowing it optionally as a synonym for "numpy" the package, but I think there is a place for "NumPy" as described above. If people like the idea, Stefan, can you put this guideline in our doc standards and examples? I'll hang it off of Developer_Zone as the start of our web style guidelines, once it's had some discussion here (and pending the outcome). --jh--

1 0