Mailman 3 March 2007 - SciPy-Dev

Unusual crash
by Russel Howe March 25, 2007

March 25, 2007

I am using today's svn scipy on OS X 10.4 intel with MacPorts python2.4. $ pythonw -c "import scipy.signal" Floating point exception $ Looking a little more closely, it is the module linalg.iterative that causes the crash. A gdb backtrace is below, has anybody else run into this before? Any ideas? I rebuilt python and scipy with no change. $ gdb /opt/local/Library/Frameworks/Python.framework/Versions/2.4/ Resources/Python.app/Contents/MacOS/Python <snip gdb banner> (gdb) run -c "import scipy.linalg.iterative" <snip loading libraries...> Program received signal EXC_ARITHMETIC, Arithmetic exception. 0x900dfb26 in strtod_l () (gdb) bt #0 0x900dfb26 in strtod_l () #1 0x90030a66 in strtod () #2 0x002ac1be in PyOS_ascii_strtod (nptr=0xbfff8680 "1.", '0' <repeats 15 times>, "1e-05", endptr=0x0) at Python/pystrtod.c:159 #3 0x002ac469 in PyOS_ascii_atof (nptr=0xbfff8680 "1.", '0' <repeats 15 times>, "1e-05") at Python/pystrtod.c:257 #4 0x002a04c7 in r_object (p=0xbfffcd50) at Python/marshal.c:484 #5 0x002a0782 in r_object (p=0xbfffcd50) at Python/marshal.c:581 #6 0x002a053d in r_object (p=0xbfffcd50) at Python/marshal.c:651 #7 0x002a0da1 in PyMarshal_ReadLastObjectFromFile (fp=0xa000be50) at Python/marshal.c:804 #8 0x0029d403 in load_source_module (name=0xbfffd6f7 "scipy.linalg.iterative", pathname=0xbfffd257 "/opt/local/lib/ python2.4/site-packages/scipy/linalg/iterative.py", fp=0xa000bdf8) at Python/import.c:723 #9 0x0029e233 in import_submodule (mod=0x121bf70, subname=0xbfffd704 "iterative", fullname=0xbfffd6f7 "scipy.linalg.iterative") at Python/ import.c:2266 #10 0x0029e465 in load_next (mod=0x121bf70, altmod=0x2e0160, p_name=0xbfffd704, buf=0xbfffd6f7 "scipy.linalg.iterative", p_buflen=0xbfffdaf8) at Python/import.c:2086 #11 0x0029e919 in PyImport_ImportModuleEx (name=0x12f72b4 "iterative", globals=0x57420, locals=0x57420, fromlist=0x1214d30) at Python/import.c:1921 #12 0x00271adf in builtin___import__ (self=0x0, args=0x50810) at Python/bltinmodule.c:45 #13 0x0020e213 in PyObject_Call (func=0xf260, arg=0x50810, kw=0x0) at Objects/abstract.c:1795 #14 0x00279a8a in PyEval_CallObjectWithKeywords (func=0xf260, arg=0x50810, kw=0x0) at Python/ceval.c:3430 #15 0x0027c63f in PyEval_EvalFrame (f=0x634f30) at Python/ceval.c:2020 #16 0x00280c66 in PyEval_EvalCodeEx (co=0x12fa720, globals=0x57420, locals=0x57420, args=0x0, argcount=0, kws=0x0, kwcount=0, defs=0x0, defcount=0, closure=0x0) at Python/ceval.c:2736 #17 0x00280ea6 in PyEval_EvalCode (co=0x12fa720, globals=0x57420, locals=0x57420) at Python/ceval.c:484 #18 0x0029ce94 in PyImport_ExecCodeModuleEx (name=0xbfffecb7 "scipy.linalg", co=0x12fa720, pathname=0xbfffdf6f "/opt/local/lib/ python2.4/site-packages/scipy/linalg/__init__.pyc") at Python/ import.c:631 #19 0x0029d24a in load_source_module (name=0xbfffecb7 "scipy.linalg", pathname=0xbfffdf6f "/opt/local/lib/python2.4/site-packages/scipy/ linalg/__init__.pyc", fp=0xa000bda0) at Python/import.c:909 #20 0x0029dd57 in load_package (name=0xbfffecb7 "scipy.linalg", pathname=0x121fac0 "\003") at Python/import.c:965 #21 0x0029e233 in import_submodule (mod=0x55d70, subname=0xbfffecbd "linalg", fullname=0xbfffecb7 "scipy.linalg") at Python/import.c:2266 #22 0x0029e465 in load_next (mod=0x55d70, altmod=0x55d70, p_name=0xbfffecbd, buf=0xbfffecb7 "scipy.linalg", p_buflen=0xbffff0b8) at Python/import.c:2086 #23 0x0029e95d in PyImport_ImportModuleEx (name=0x506c4 "scipy.linalg.iterative", globals=0x20a50, locals=0x20a50, fromlist=0x2e0160) at Python/import.c:1928 #24 0x00271adf in builtin___import__ (self=0x0, args=0x1fae0) at Python/bltinmodule.c:45 #25 0x0020e213 in PyObject_Call (func=0xf260, arg=0x1fae0, kw=0x0) at Objects/abstract.c:1795 #26 0x00279a8a in PyEval_CallObjectWithKeywords (func=0xf260, arg=0x1fae0, kw=0x0) at Python/ceval.c:3430 #27 0x0027c63f in PyEval_EvalFrame (f=0x608890) at Python/ceval.c:2020 #28 0x00280c66 in PyEval_EvalCodeEx (co=0x4ee20, globals=0x20a50, locals=0x20a50, args=0x0, argcount=0, kws=0x0, kwcount=0, defs=0x0, defcount=0, closure=0x0) at Python/ceval.c:2736 #29 0x00280ea6 in PyEval_EvalCode (co=0x4ee20, globals=0x20a50, locals=0x20a50) at Python/ceval.c:484 #30 0x002a80a0 in PyRun_SimpleStringFlags (command=0x600320 "import scipy.linalg.iterative\n", flags=0xbffff618) at Python/pythonrun.c:1265 #31 0x002b1493 in Py_Main (argc=1, argv=0xbffff6a0) at Modules/main.c: 481 #32 0x000018ee in ?? () #33 0x00001815 in ?? () (gdb)

2 7

New Operators in Python
by Travis Oliphant March 24, 2007

March 24, 2007

Every so often the idea of new operators comes up because of the need to do both "matrix-multiplication" and element-by-element multiplication. I think this is one area where the current Python approach is not as nice because we have a limited set of operators to work with. One thing I wonder is if we are being vocal enough with the Python 3000 crowd to try and get additional operators into the language itself. What if we could get a few new operators into the language to help us. If we don't ask for it, it certainly won't happen. My experience is that the difficulty of using the '*' operator for both matrix multiplication and element-by-element multiplication depending on the class of the object is not especially robust. It makes it harder to write generic code, and we still haven't gotten everything completely right. It is somewhat workable as it stands, but I think it would be nicer if we could have some "meta" operator that allowed an alternative definition of major operators. Something like @* for example (just picking a character that is already used for decorators). I wonder if we should propose such a thing for Python 3000. -Travis

7 11

Have divide by zero handling changed lately?
by Otto Tronarp March 22, 2007

March 22, 2007

Hello, have numpy's handling of divide by zero changed lately? import numpy numpy.seterr(divide='warn') v = numpy.asarray([0, 1], dtype='d') print 1/numpy.asarray([0], dtype='d') print 1/v print 1.0/v[0] In the code above 1.0/v[0] doesn't generate a warning, that can't be right? While writing this I actually tried some more things: 1.0/numpy.float32(0) has correct behavior (generates a warning) whereas 1.0/numpy.float64(0) silently gives us an inf. This is with In [95]:numpy.__version__ Out[95]:'1.0.1' on a i686 AMD Sempron(tm) 2600+ AuthenticAMD GNU/Linux Regards, Otto

1 0

Re: [SciPy-dev] linsolve.factorized was: Re: Using umfpack to calculate a incomplete LU factorisation (ILU)
by Neilen Marais March 22, 2007

March 22, 2007

Hi Robert! On Thu, 08 Mar 2007 18:10:30 +0100, Robert Cimrman wrote: > Robert Cimrman wrote: > Well, I did it since I am going to need this, too :-) > > In [3]:scipy.linsolve.factorized? > ... > Definition: scipy.linsolve.factorized(A) > Docstring: > Return a fuction for solving a linear system, with A pre-factorized. > > Example: > solve = factorized( A ) # Makes LU decomposition. > x1 = solve( rhs1 ) # Uses the LU factors. > x2 = solve( rhs2 ) # Uses again the LU factors. > > This uses UMFPACK if available. This is a useful improvement, thanks. But why not just extend linsolve.splu to use umfpack so we can present a consistent interface? The essential difference between factorized and splu is that you get to explicity control the storage of the LU factorisation and get some additional info (i.e. the number of nonzeros), whereas factorised only gives you a solve function. The actual library used to do the sparse LU is just an implementation detail that should abstracted wherever possible, no? If nobody complains about the idea I'm willing to implement it. Thanks Neilen > > cheers, > r. -- you know its kind of tragic we live in the new world but we've lost the magic -- Battery 9 (www.battery9.co.za)

2 1

[Fwd: Re: wrapping BLOPEX for SciPy]
by Robert Cimrman March 22, 2007

March 22, 2007

Hi all, FYI: I have asked Andrew Knyazev, the author of BLOPEX (sparse) eigenvalue solver, some questions related to making it usable from SciPy (which I plan to do, unless somebody else did/does it already). Some interesting points have appeared in the discussion, see below remarks of Julien Langou, one of the LAPACK developers (the discussion leading to those remarks is below below :). r. -------- Original Message -------- Subject: Re: wrapping BLOPEX for SciPy Date: Wed, 21 Mar 2007 14:28:24 -0600 (MDT) From: Julien Langou <langou(a)math.cudenver.edu> To: Andrew Knyazev <andrew.knyazev(a)cudenver.edu> CC: Robert Cimrman <cimrman3(a)ntc.zcu.cz> Hello, Some points from the top of my heads for some SciPY developers. * while Matlab eig correctly call SYEV (symmetric eigenvalue solver) or GEEV (nonsymmetric eigenvalue solver) from LAPACK depending whether the matrix is symmetric or not, Matlab eigs call ARPACK whether the matrix is symmetric or not. I always thought this was a big liability for such a software. Symmetric eigensolver would be way faster, more stable, eigenvectors would have guaranteed orthogonality and eigenvalues would be real. Below a Matlab script that shows that when there is a cluster of eigenvalue, eigs (ARPACK) returns the correct invariant subspace but the eigenvectors are not orthogonalized. (logic for a nonsymmetric eigensolver ...) >> clear >> n=10; V=orth(randn(n)); D=1:n; D(n-1)=n; A=V*diag(D)*V'; >> opts.disp=0; >> [Veigs,Deigs] = eigs(A,2,'LM',opts); >> Deigs Deigs = 10 0 0 10 >> Veigs'*Veigs ans = 1.0000 0.9545 0.9545 1.0000 Conclusion having a true symmetric eigensolver (e.g. lobpcg) for the SciPy eigs would be clearly a big plus with respect to Matlab. * They are four LAPACK symmetric eigensolvers. SYEV (symmetric QR), SYEVD (divide and conquer), SYEVD (divide and conquer), SYEVX (bisection and inverse iteration) and SYEVR (MRRR). Matlab uses SYEV. The most robust but by far the slowest. We have worked on all those algoirthms and they are now all robust. For MRRR, this is new. The version in 3.1 is really stable, the one in 3.0, well .... In term of performance, there is no comparison SYEVR and SYEVX are way faster than SYEV. Performance of SYEVX are really matrix dependent. Comparing SYEVD and SYEVR is hard. * SYEVX and SYEVR supports subset computation (same thing: subset computation in SYEVR is from 3.1). So in some sense if you really want to compute only 10 of the eigenvalue of a matrix with LAPACK, you can save a good bunch of FLOPS. This option is not available from Matlab. Note that the cost is still O(n^3) (due to tridiagonalization). * LAPACK-3.1 nonsymmetric eigensolver is about 3-5 times faster than the LAPACK-3.0 one. Matlab has still not yet upgraded to 3.1, as far as I know. * For the question whether to use eigs or eig; or equivallently when to use ARPACK or LAPACK. This depend on the size of the matrix, the sparsity of the matrix and the number of eigenvalues desired. Just to check attached is a plot of comparison betwenn eig and eigs for dense matrices going from n=100:100:1000, for k=[1 10]. Conclusion: eigs behaves as an N^2 algorithm as expected, eig behaves as a N^3 algorithm as expected. On my machine, if you want one eigenvalue, you better off with LAPACK for n<150. If you want 10 eigenvalues, you better off with LAPACK for n<450. Note that the scale is logarithmic so when n=1000, the n^3 factor completely kills LAPACK and it is about 5 times slower than eigs. -j On Wed, 21 Mar 2007, Andrew Knyazev wrote: > Robert Cimrman wrote: >> Andrew Knyazev wrote: >>> LAPACK is indeed for full matrices only. If you want all >>> eigenvectors, the matrix of all eigenvectors is usually full even >>> if your original matrix is sparse. Since you loose the sparsity >>> at the end anyway, the standard procedure is to convert your >>> sparse input matrix into full and plug it into LAPACK. Or may be >>> you try to do something non orthodox that I do not know. >> >> OK. What SciPy needs is an equivalent of Matlab's eigs. I have >> briefly looked at it, and it switches to eig( full( A ) ) only if >> all eigenvalues are requested. I would say that if one wants e.g. k >> = 999 eigenvalues out of n = 1000, making a full matrix would still >> be preferable. Is there any recommendation on for which k to switch >> to the full matrix mode, especially when using BLOPEX? >> >> Thank you for your time and answers! >> >> r. > > I never did a comparison between MATLAB's eig (which is an interface > to LAPACK) and eigs (which is an interface to ARPACK). I would guess > that under normal circumstances, e.g., having enough memory, eig > should be faster than eigs if >50% eigenvectors are requested, may be > less. With BLOPEX, I would guess that eig is faster for >10% (if no > preconditioning is used). The main advantage of BLOPEX is that it > allows using preconditioning to accelerate the convergence. The > preconditioning must be provided by the user, though. > > I copy Julien Langou, one of the LAPACK developers, who may be able > to add some more to this. >

1 0

confidence intervals on multi-parameter minimisation
by Graeme O'Keefe March 21, 2007

March 21, 2007

I have been using fmin_l_bfgs_b to perform bounded non-linear multi- parameter estimation. Is there a way to determine the confidence intervals on the resultant parameter estimates. I've seen reference to scipy.odr as a better alternative to scipy.optimize and have checked out the latest svn and currently trying to build (OS X 10.4.9) which is always fun. However, I can't find any documentation on what scipy.odr is. Can someone point me somewhere. regards, Graeme

3 7

Google Summer of Code 2007
by Jarrod Millman March 20, 2007

March 20, 2007

Hey, I was wondering if Enthought would be interested in signing up as a mentor organization: http://groups.google.com/group/google-summer-of-code-announce/web/gsoc-ment… It is a fair amount of work. Fernando Perez explained that PSF wants to focus on core Python development this year. More specifically after google tells them how many slots they get, the PSF plans to reserve a certain number of slots set aside specifically for core projects. Given the increase in interest in the summer of code, the PSF may not have many slots for scientific python. The reason this is important is that I have been talking with two students, David Cournapeau and Taylor Berg, about potential projects. They would both be related to machine learning and I think would be very useful contributions to SciPy. Fernando Perez and Brian Granger are also thinking about mentoring Min on an IPython project. I don't know if anyone else is thinking about other possible SciPy-related SoC projects. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

4 5

New Trac feature: TracReSTMacro
by Jeff Strunk March 19, 2007

March 19, 2007

Good afternoon, By request, I have installed the TracReSTMacro on the numpy, scipy, and scikits tracs. This plugin allows you to display ReST formatted text directly from svn. For example, http://projects.scipy.org/neuroimaging/ni/wiki/ReadMe in its entirety is: [[ReST(/ni/trunk/README)]] Thank you, Jeff Strunk IT Administrator Enthought, Inc.

1 0

inconsistent use of tabs and spaces in indentation in numpy\distutils\log.py
by Berthold Höllmann March 19, 2007

March 19, 2007

Calling my setup script with the -tt option to python I get: Traceback (most recent call last): File "setup.py", line 24, in ? from numpy.distutils.core import setup File "C:\Python24\lib\site-packages\numpy\distutils\__init__.py", line 5, in ? import ccompiler File "C:\Python24\lib\site-packages\numpy\distutils\ccompiler.py", line 11, in ? import log File "C:\Python24\lib\site-packages\numpy\distutils\log.py", line 7, in ? from misc_util import red_text, yellow_text, cyan_text, is_sequence, is_string File "C:\Python24\lib\site-packages\numpy\distutils\misc_util.py", line 975 import numpy.numarray.util as nnu ^ TabError: inconsistent use of tabs and spaces in indentation I would report this as bug, but numpy.scipy.org shows only the fedora apache default screen. This is numpy 1.0.1 Kind regards Berthold Höllmann -- Germanischer Lloyd AG CAE Development Vorsetzen 35 20459 Hamburg Phone: +49(0)40 36149-7374 Fax: +49(0)40 36149-7320 e-mail: berthold.hoellmann(a)gl-group.com Internet: http://www.gl-group.com This e-mail and any attachment thereto may contain confidential information and/or information protected by intellectual property rights for the exclusive attention of the intended addressees named above. Any access of third parties to this e-mail is unauthorised. Any use of this e-mail by unintended recipients such as total or partial copying, distribution, disclosure etc. is prohibited and may be unlawful. When addressed to our clients the content of this e-mail is subject to the General Terms and Conditions of GL's Group of Companies applicable at the date of this e-mail. If you have received this e-mail in error, please notify the sender either by telephone or by e-mail and delete the material from any computer. GL's Group of Companies does not warrant and/or guarantee that this message at the moment of receipt is authentic, correct and its communication free of errors, interruption etc. Germanischer Lloyd AG, 31393 AG HH, Hamburg, Vorstand: Dr. Hermann J. Klein, Rainer Schöndube, Vorsitzender des Aufsichtsrats: Dr. Wolfgang Peiner

2 3

ERROR: test_explicit (scipy.tests.test_odr.test_odr)
by Nils Wagner March 15, 2007

March 15, 2007

Hi all, scipy.test(1) results in 5 errors wrt odr package. It seems to be a python2.5 problem since I cannot reproduce these errors with python2.4 and >>> scipy.__version__ '0.5.3.dev2848' Anyway can someone reproduce these errors ? ERROR: test_explicit (scipy.tests.test_odr.test_odr) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/odr/tests/test_odr.py", line 46, in test_explicit out = explicit_odr.run() File "/usr/local/lib64/python2.5/site-packages/scipy/odr/odrpack.py", line 1049, in run self.output = Output(apply(odr, args, kwds)) File "/usr/local/lib64/python2.5/site-packages/scipy/odr/odrpack.py", line 576, in __init__ self.stopreason = report_error(self.info) File "/usr/local/lib64/python2.5/site-packages/scipy/odr/odrpack.py", line 143, in report_error 'Iteration limit reached')[info % 10] IndexError: tuple index out of range Nils

1 0