Mailman 3 September 2012 - SciPy-User

Avoiding inner for loops??
by Martin De Kauwe Dec. 22, 2012

Dec. 22, 2012

Hi, I need to avoid (at least) two inner for loops in what I am trying to do otherwise my processing takes forever. What is the best way to transfer what I am doing into a more "numpy way"? Essentially I am trying to call a model again for various different parameter combinations. The example is fictional, by the grid_size would ideally grow > 500 and by doing so the processing speed becomes very slow the way I have set things up.. thanks. example. import numpy as np def fake_model(data1, data2, p1, p2, p3): """ complete nonsense """ return data1 + data2 * p1 * p2 * p3 data1 = np.random.rand(10) # the size of this arrays varies might be 10 might be 15 etc data2 = np.random.rand(10) # the size of this arrays varies might be 10 might be 15 etc obs = np.random.rand(10) # the size of this arrays varies might be 10 might be 15 etc grid_size = 10 # Ideally this would be a large number param1 = np.linspace(5.0, 350, grid_size) param2 = np.linspace(5.0, 550, grid_size) param3 = np.linspace(1E-8, 10.5, grid_size) ss = np.zeros(0) for p1 in param1: for p2 in param2: for p3 in param3: ans = fake_model(data1, data2, p1, p2, p3) ss = np.append(ss, np.sum(obs - ans)**2) print np.sum(obs - ans)**2

6 12

Request for usage examples on scipy.stats.rv_continuous and scipy.ndimage.grey_dilate(structure)
by Christoph Deil Dec. 5, 2012

Dec. 5, 2012

Hi, I am having some trouble figuring out how to use the following two aspects of scipy: 1) How can I create a random distribution given by some python formula, e.g. pdf(x) = x**2 in the range x = 1 to 2. The example given in the rv_continuous docstring doesn't work for me: In [1]: import matplotlib.pyplot as plt In [2]: numargs = generic.numargs AttributeError: type object 'numpy.generic' has no attribute 'numargs' Is it also possible to use rv_continuous for multidimensional distributions, say pdf(x,y)=x*y in the range x = 1 to 2 and y = 0 to 3? Or for pdfs that are just given numerically on a grid (say as the result of a gaussian_kde() ), not by a python function? It would be nice if someone could add a few examples that show rv_continuous usage to the docstring or tutorial. 2) How to use the 'structure' argument of scipy.ndimage.grey_dilate() correctly? There are some explanations at http://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html but I think they are inconsistent with the implementation: -------------------------------------------------------------- In [98]: a = zeros((7,7)); a[3,3]=1; a Out[98]: array([[ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 1., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.]]) In [100]: s = generate_binary_structure(2,2); s Out[100]: array([[ True, True, True], [ True, True, True], [ True, True, True]], dtype=bool) In [101]: grey_dilation(a,size=(3,3)) Out[101]: array([[ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.]]) In [102]: grey_dilation(a,size=(3,3),structure=s) Out[102]: array([[ 1., 1., 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1., 1., 1.], [ 1., 1., 2., 2., 2., 1., 1.], [ 1., 1., 2., 2., 2., 1., 1.], [ 1., 1., 2., 2., 2., 1., 1.], [ 1., 1., 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1., 1., 1.]]) In [103]: grey_dilation(a,structure=s) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /Users/deil/<ipython console> in <module>() /opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/ndimage/morphology.pyc in grey_dilation(input, size, footprint, structure, output, mode, cval, origin) 354 sz = footprint.shape[ii] 355 else: --> 356 sz = size[ii] 357 if not sz & 1: 358 origin[ii] -= 1 TypeError: 'NoneType' object is unsubscriptable In [104]: grey_dilation(a,footprint=s) Out[104]: array([[ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 1., 1., 1., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0., 0., 0.]]) -------------------------------------------------------------- Why doesn't grey_dilation(a,structure=s) work? Why doesn't grey_dilation(a,size=(3,3),structure=s) give the same result as grey_dilation(a,size=(3,3),footprint=s)? I am running on numpy version 1.4.0 and scipy version 0.7.1 Thanks! --Christoph

8 11

Problem with Installation of Scipy on Macbook
by manisha pujari Nov. 21, 2012

Nov. 21, 2012

Hello everyone, This is the first time I am using scipy and I am having much trouble to just install it on my macbook. I am using Python 2.6 and numpy 1.6.0. I downloaded scipy0.9.0 tar.gz file from the web link http://sourceforge.net/projects/scipy/files/scipy/ and tried to install it from the source. But on giving the command scipy-0.9.0$ python setup.py build it is always giving the following problem : Warning: No configuration returned, assuming unavailable.blas_opt_info: FOUND: extra_link_args = ['-Wl,-framework', '-Wl,Accelerate'] define_macros = [('NO_ATLAS_INFO', 3)] extra_compile_args = ['-faltivec', '-I/System/Library/Frameworks/vecLib.framework/Headers'] non-existing path in 'scipy/io': 'docs' lapack_opt_info: FOUND: extra_link_args = ['-Wl,-framework', '-Wl,Accelerate'] define_macros = [('NO_ATLAS_INFO', 3)] extra_compile_args = ['-faltivec'] umfpack_info: libraries umfpack not found in /System/Library/Frameworks/Python.framework/Versions/2.6/lib libraries umfpack not found in /usr/local/lib libraries umfpack not found in /usr/lib /Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/system_info.py:463: UserWarning: UMFPACK sparse solver (http://www.cise.ufl.edu/research/sparse/umfpack/) not found. Directories to search for the libraries can be specified in the numpy/distutils/site.cfg file (section [umfpack]) or by setting the UMFPACK environment variable. warnings.warn(self.notfounderror.__doc__) NOT AVAILABLE Traceback (most recent call last): File "setup.py", line 181, in <module> setup_package() File "setup.py", line 173, in setup_package configuration=configuration ) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/core.py", line 152, in setup config = configuration() File "setup.py", line 122, in configuration config.add_subpackage('scipy') File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 972, in add_subpackage caller_level = 2) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 941, in get_subpackage caller_level = caller_level + 1) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 878, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/setup.py", line 20, in configuration config.add_subpackage('special') File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 972, in add_subpackage caller_level = 2) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 941, in get_subpackage caller_level = caller_level + 1) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 878, in _get_configuration_from_setup_py config = setup_module.configuration(*args) File "scipy/special/setup.py", line 54, in configuration extra_info=get_info("npymath") File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 2184, in get_info pkg_info = get_pkg_info(pkgname, dirs) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/misc_util.py", line 2136, in get_pkg_info return read_config(pkgname, dirs) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/npy_pkg_config.py", line 390, in read_config v = _read_config_imp(pkg_to_filename(pkgname), dirs) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/npy_pkg_config.py", line 326, in _read_config_imp meta, vars, sections, reqs = _read_config(filenames) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/npy_pkg_config.py", line 309, in _read_config meta, vars, sections, reqs = parse_config(f, dirs) File "/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/distutils/npy_pkg_config.py", line 281, in parse_config raise PkgNotFound("Could not find file(s) %s" % str(filenames)) numpy.distutils.npy_pkg_config.PkgNotFound: Could not find file(s) ['/Library/Python/2.6/site-packages/numpy-1.6.0-py2.6-macosx-10.6-universal.egg/numpy/core/lib/npy-pkg-config/npymath.ini'] I am unable to understand where and what exactly the problem is. Can anyone please help me? I will be really thankful for a response. -- Regards, Manisha

4 3

ANN: WinPython v2.7.3.0
by Pierre Raybaut Oct. 11, 2012

Oct. 11, 2012

Hi all, I'm pleased to introduce my new contribution to the Python community: WinPython. WinPython v2.7.3.0 has been released and is available for 32-bit and 64-bit Windows platforms: http://code.google.com/p/winpython/ WinPython is a free open-source portable distribution of Python for Windows, designed for scientists. It is a full-featured (see http://code.google.com/p/winpython/wiki/PackageIndex) Python-based scientific environment: * Designed for scientists (thanks to the integrated libraries NumPy, SciPy, Matplotlib, guiqwt, etc.: * Regular *scientific users*: interactive data processing and visualization using Python with Spyder * *Advanced scientific users and software developers*: Python applications development with Spyder, version control with Mercurial and other development tools (like gettext) * *Portable*: preconfigured, it should run out of the box on any machine under Windows (without any installation requirements) and the folder containing WinPython can be moved to any location (local, network or removable drive) * *Flexible*: one can install (or should I write "use" as it's portable) as many WinPython versions as necessary (like isolated and self-consistent environments), even if those versions are running different versions of Python (2.7, 3.x in the near future) or different architectures (32bit or 64bit) on the same machine * *Customizable*: using the integrated package manager (wppm, as WinPython Package Manager), it's possible to install, uninstall or upgrade Python packages (see http://code.google.com/p/winpython/wiki/WPPM for more details on supported package formats). *WinPython is not an attempt to replace Python(x,y)*, this is just something different (see http://code.google.com/p/winpython/wiki/Roadmap): more flexible, easier to maintain, movable and less invasive for the OS, but certainly less user-friendly, with less packages/contents and without any integration to Windows explorer [*]. [*] Actually there is an optional integration into Windows explorer, providing the same features as the official Python installer regarding file associations and context menu entry (this option may be activated through the WinPython Control Panel). Enjoy! -Pierre

3 2

Bessel function of complex order
by Andreas Pritschet Oct. 5, 2012

Oct. 5, 2012

Hi, I have noticed in the docs and some "bug reports" that Bessel functions in SciPy support only real order. But for my work I require a modified Bessel function of second kind of complex(!) order for complex values. Is in SciPy a chance of calculating something like scipy.special.kv(1j*k,1j), whereby k is an array?? Thanks and best regards Andi -- Andreas Pritschet Phone: +49 151 11728439 Homepage: http://www.pritschet.me GPG Pub Key: http://goo.gl/4mOsM

5 4

Fitting Gaussian in spectra
by Joe Philip Ninan Oct. 4, 2012

Oct. 4, 2012

Hi, I have a spectra with multiple gaussian emission lines over a noisy continuum. My primary objective is to find areas under all the gaussian peaks. For that, the following is the algorithm i have in mind. 1) fit the continuum and subtract it. 2) find the peaks 3) do least square fit of gaussian at the peaks to find the area under each gaussian peaks. I am basically stuck at the first step itself. Simple 2nd or 3rd order polynomial fit is not working because the contribution from peaks are significant. Any tool exist to fit continuum ignoring the peaks? For finding peaks, i tried find_peaks_cwt in signal module of scipy. But it seems to be quite sensitive of the width of peak and was picking up non-existing peaks also. The wavelet used was default mexican hat. Is there any better wavelet i should try? Or is there any other module in python/scipy which i should give a try? Thanking you. -cheers joe -- /--------------------------------------------------------------- "GNU/Linux: because a PC is a terrible thing to waste" - GNU Generation

10 10

Re: [SciPy-User] Pylab - standard packages
by William Furnass Oct. 2, 2012

Oct. 2, 2012

Something that could raise awareness and expedite the adoption of the Pylab standard under Linux would be the availability of repositories for some of the most common distributions (e.g. Ubuntu, Fedora, RHEL-a-likes). As a number of distributions take a while to catch up with the latest releases of ipython etc making Ubuntu PPA and repos.fedorapeople.org repositories available could provide people with a familar, quick and easy means to install Pylab. Not sure whether we would want distribution package managers automatically upgrading packages though when new Pylab standards/packages are released. Any thoughts? On a related note it might be a good opportunity to bring the fairly-official-looking Scipy PPA [1] up to date. [1] https://launchpad.net/~scipy/+archive/ppa Cheers, Will

9 15

Request help on fsolve
by The Helmbolds Oct. 1, 2012

Oct. 1, 2012

Please help me out here. I’m trying to rewrite the docstring for the `fsolve.py` routine located on my machine in: C:/users/owner/scipy/scipy/optimize/minpack.py The specific issue I’m having difficulty with is understanding the outputs described in fsolve’s docstring as: 'fjac': the orthogonal matrix, q, produced by the QR factorization of the final approximate Jacobian matrix, stored column wise 'r': upper triangular matrix produced by QR factorization of same matrix These are described in SciPy’s minpack/hybrd.f file as: ‘fjac’ is an output n by n array which contains the orthogonal matrix q produced by the qr factorization of the final approximate jacobian. ‘r’ is an output array of length lr which contains the upper triangular matrix produced by the qr factorization of the final approximate jacobian, stored rowwise. For ease in writing, in what follows let’s use the symbols ‘Jend’ for the final approximate Jacobian matrix, and use ‘Q’ and ‘R’ for its QR decomposition matrices. Now consider the problem of finding the solution to the following three nonlinear equations in three unknowns (u, v, w), which we will refer to as ‘E’: 2 * a * u + b * v + d - w * v = 0 b * u + 2 * c * v + e - w * u = 0 u * v - f = 0 where (a, b, c, d, e, f ) = (2, 3, 7, 8, 9, 2). For inputs to fsolve, we identify (u, v, w) = (x[0], x[1], x[2]). Now fsolve gives the solution array: [uend vend wend] = [ 1.79838825 1.11210691 16.66195357]. With these values, the above three equations E are satisfied to an accuracy of about 9 significant figures. The Jacobian matrix for the three LHS functions in E is: J = np.matrix([[2*a, b-w, -v], [b-w, 2*c, -u], [v, u, 0.]]) Note that it’s symmetric, and if we compute its value using the above fsolve’s ‘end’ solution values we get: Jend = [[ 4. 19.66195357 1.11210691], [ 19.66195357 14. 1.79838825], [ 1.11210691 1.79838825 0. ]] Using SciPy’s linalg package, this Jend has the QR decomposition: Qend = [[-0.28013447 -0.91516674 -0.28981807] [ 0.95679602 -0.24168763 -0.16164302] [ 0.07788487 -0.32257856 0.94333293]] Rend = [[-14.278857 17.08226116 -1.40915124] [ -0. 9.69946027 1.45241144] [ -0. 0. 0.61300558]] and Qend * Rend = Jend to within about 15 significant figures. However, fsolve gives the QR decomposition: qretm = [[-0.64093238 0.75748326 0.1241966 ] [-0.62403598 -0.60841098 0.4903215 ] [-0.44697291 -0.23675978 -0.8626471 ]] rret = [ -7.77806716 30.02199802 -0.819055 -10.74878184 2.00090268 1.02706198] and converting rret to a NumPy matrix gives: rretm = [[ -7.77806716 30.02199802 -0.819055 ] [ 0. -10.74878184 2.00090268] [ 0. 0. 1.02706198]] Now qret and rretm bear no obvious relation to Qend and Rend. Although qretm is orthogonal to about 16 significant figures, we find the product: qretm * rretm = [[ 4.98521509 -27.38409295 2.16816676] [ 4.85379376 -12.19513008 -0.2026608 ] [ 3.47658529 -10.87414051 -0.99362993]] which bears no obvious relationship to Jend. The hybrdj.f routine in minpack refers to a permutation matrix, p, such that we should have in our notation: p*Jend = qretm*rretm, but fsolve apparently does not return the matrix p, and I don’t see any permutation of Jend that would equal qretm*rretm. If we reinterpret rret as meaning the matrix: rretaltm = [[ -7.77806716 30.02199802 -10.74878184] [ 0. -0.819055 2.00090268] [ 0. 0. 1.02706198]] then we get the product: qretm * rretaltm = [[ 4.98521509 -19.86249109 8.53245022] [ 4.85379376 -18.2364849 5.99384603] [ 3.47658529 -13.22510045 3.44468895]] which again bears no obvious relationship to Jend. Using the transpose of qretm in the above product is no help. So please help me out here. What are the fjac and r values that fsolve returns? How are they related to the above Qend, Rend, and Jend? How is the user supposed to use them? Bob H

2 1

ANN: scikits-image 0.7.0 release
by Tony Yu Sept. 30, 2012

Sept. 30, 2012

Announcement: scikits-image 0.7.0 ================================= We're happy to announce the 7th version of scikits-image! Scikits-image is an image processing toolbox for SciPy that includes algorithms for segmentation, geometric transformations, color space manipulation, analysis, filtering, morphology, feature detection, and more. For more information, examples, and documentation, please visit our website http://skimage.org New Features ------------ It's been only 3 months since scikits-image 0.6 was released, but in that short time, we've managed to add plenty of new features and enhancements, including - Geometric image transforms - 3 new image segmentation routines (Felsenzwalb, Quickshift, SLIC) - Local binary patterns for texture characterization - Morphological reconstruction - Polygon approximation - CIE Lab color space conversion - Image pyramids - Multispectral support in random walker segmentation - Slicing, concatenation, and natural sorting of image collections - Perimeter and coordinates measurements in regionprops - An extensible image viewer based on Qt and Matplotlib, with plugins for edge detection, line-profiling, and viewing image collections Plus, this release adds a number of bug fixes, new examples, and performance enhancements. Contributors to this release ---------------------------- This release was only possible due to the efforts of many contributors, both new and old. - Andreas Mueller - Andreas Wuerl - Andy Wilson - Brian Holt - Christoph Gohlke - Dharhas Pothina - Emmanuelle Gouillart - Guillaume Gay - Josh Warner - James Bergstra - Johannes Schonberger - Jonathan J. Helmus - Juan Nunez-Iglesias - Leon Tietz - Marianne Corvellec - Matt McCormick - Neil Yager - Nicolas Pinto - Nicolas Poilvert - Pavel Campr - Petter Strandmark - Stefan van der Walt - Tim Sheerman-Chase - Tomas Kazmar - Tony S Yu - Wei Li

1 0

Re: [SciPy-User] COBYLA
by The Helmbolds Sept. 30, 2012

Sept. 30, 2012

Paul Virtanen writes: >>Date: Thu, 27 Sep 2012 23:57:42 +0300 >>From: Pauli Virtanen <pav(a)iki.fi> >>Subject: Re: [SciPy-User] Flaws in fmin_cobyla (SNIP) >>Ideally, Ipython consoles et al. would redirect stdout/stderr file >>handles so that also the output that comes outside from Python land >>would be shown to the user. However, the correct behavior for Scipy >>routines would be to print via Python's I/O. That's OK for "one-shot" optimizations. Not so good for real-world "mutli-shot" optimizations -- unless I'm missing something here. By multi-shot optimizations I mean something like the following: An engine designer has an 18-dimensional design variable, with 24 constraints and 60 args, and wants to run 70 to 100 cases by varying the args. She writes a Python program that reads the args from a file of cases, prepares the inputs to cobyla, turns the cobyla crank and saves its outputs to a file for further analysis. She wants to save all of the Results information from each run to this output file. That's not convenient to do if the Results information is not actually _returned_ by COBYLA. Bob and Paula H "Many of life's failures are people who did not realize how close they were to success when they gave up." (Thomas Edison)

3 2