Mailman 3 November 2006 - SciPy-Dev

ARPACK wrapper
by Neilen Marais Nov. 21, 2006

Nov. 21, 2006

Hi I've commited the beginnings of an ARPACK wrapper partially adressing http://projects.scipy.org/scipy/scipy/ticket/231 in r2323 I based it on the arpack-0.10.tar.bz2 source that aric attached to #231. I added an ARPACK driver that demonstrates how to use ARPACK's shift-invert mode to solve a generalised eigensystem. Unfortunately I haven't generated any compact test data for that yet (I was using matrices created by my own FEM app) but will soon. As things stand now a lot of work still remains to be done. I can think of (off the top of my head) 1) Merging Aric and my drivers 2) Handling more cases with an easy-to-use interface 3) Using the symmetric/complex/etc. modes 4) Many more, no doubt ;) I'm away for the rest of the weekend, but I look forward to comments on Monday :) Regards Neilen -- you know its kind of tragic we live in the new world but we've lost the magic -- Battery 9 (www.battery9.co.za)

3 8

ARPACK Wrapper update
by Neilen Marais Nov. 20, 2006

Nov. 20, 2006

Hi. I spent some time looking at Nils's problem (I'll follow up on his post), and also refactoring the wrappers a bit. First off, I'd like to get some opinions. I think we should structure the modules such that a user doesn't know they are dealing with ARPACK, but merely a sparse eigensolver unless they want to do advanced things. This should be in a module called speigs. It should present an interface such as Aric's code in arpack.py. Then there should be a medium-level interface that takes care of ARPACK's Fortran nastiness for you but let's you use more advanced features. This is what is currently in the file called speigs.py. The user-level code should call the medium-level ARPACK interface to get the work done. And of course the raw Fortran wrappers should be available to control freaks. The way stuff is now (which is purely a coincidence and not really by design) is backwards. My code that is now in speigs.py should be in arpack.py, and vice-versa. I refactored the code in speigs.py so that the generalised and ordinary solvers don't duplicade code. I also renamed the function calls to start with ARPACK to make it clear that they are tied to ARPACK and aren't general user routines. So, I propose to: 1) rename speigs.py to arpack.py, and vice-versa 2) Make speigs and arpack separate modules 3) Modify the code in the proposed speigs.py to call the code in the proposed arpack.py to do the actual ARPACK work Does this sound good? At the moment the ARPACK_ functions can handle real, general, double-precision matrices, and use the following ARPACK modes: a) Generalised eigen-problems with spectrum shift b) Standard eigen-problems with no spectrum shift It should also support: c) Generalised eigen-problems with no spectrum shift d) Standard with spectrum shift Adding c) and d) should not be too much work. Once that is done we should look for a reasonable way to support complex matrices and other numerical types. I think Aric's wrappers already address that to some extent, so I'll take a look there. Thoughts? Regards Neilen -- you know its kind of tragic we live in the new world but we've lost the magic -- Battery 9 (www.battery9.co.za)

1 0

Current status of arpack.speigs.eigvals
by Nils Wagner Nov. 20, 2006

Nov. 20, 2006

Hi all, I am wondering if arpack.speigs.eigvals is ready to compute eigenvalues. a = random.rand(10,10) matvec = lambda x:dot(a,x) nev = 4 ncv = 2*nev w1,v1 = arpack.speigs.eigvals(matvec,a.shape[0],nev=nev,ncv=ncv) Traceback (most recent call last): File "test_eigs.py", line 34, in ? w1,v1 = arpack.speigs.eigvals(matvec,a.shape[0],nev=nev,ncv=ncv) File "/usr/lib64/python2.4/site-packages/scipy/sandbox/arpack/speigs.py", line 71, in eigvals if info != 0: raise "Hell" # Indicates some error during the Arnouldi iterations Hell Am I missing something ? Nils

2 1

Re: [SciPy-dev] PyEM problem
by David Cournapeau Nov. 16, 2006

Nov. 16, 2006

Axel Brink wrote: > Dear David, > > I run into a problem when trying your PyEM example "Creating, sampling > and plotting a mixture": > > Python 2.4.3 (#1, Jun 13 2006, 11:46:08) > [GCC 4.1.1 20060525 (Red Hat 4.1.1-1)] on linux2 > Type "help", "copyright", "credits" or "license" for more information. > >>> import numpy as N > >>> import pylab as P > >>> from scipy.sandbox.pyem import GM > using scipy.cluster.vq > >>> > >>> #------------------------------ > ... # Hyper parameters: > ... # - K: number of clusters > ... # - d: dimension > ... k = 3 > >>> d = 2 > >>> > >>> #------------------------------------------------------- > ... # Values for weights, mean and (diagonal) variances > ... # - the weights are an array of rank 1 > ... # - mean is expected to be rank 2 with one row for one component > ... # - variances are also expteced to be rank 2. For diagonal, one row > ... # is one diagonal, for full, the first d rows are the first > variance, > ... # etc... In this case, the variance matrix should be k*d rows and d > ... # colums > ... w = N.array([0.2, 0.45, 0.35]) > >>> mu = N.array([[4.1, 3], [1, 5], [-2, -3]]) > >>> va = N.array([[1, 1.5], [3, 4], [2, 3.5]]) > >>> > >>> #----------------------------------------- > ... # First method: directly from parameters: > ... # Both methods are equivalents. > ... gm = GM.fromvalues(w, mu, va) > >>> > >>> #------------------------------------- > ... # Second method to build a GM instance: > ... gm = GM(k, d, mode = 'diag') > >>> # The set_params checks that w, mu, and va corresponds to k, d and m > ... gm.set_params(w, mu, va) > Traceback (most recent call last): > File "<stdin>", line 2, in ? > AttributeError: GM instance has no attribute 'set_params' > > > I also tried 'set_param', but this doesn't help either: > > >>> gm.set_param(w, mu, va) > Traceback (most recent call last): > File "<stdin>", line 1, in ? > File > "/usr/lib/python2.4/site-packages/scipy/sandbox/pyem/gauss_mix.py", > line 84, in set_param > raise GmParamError("Number of given components is %d, expected %d" > NameError: global name 'shape' is not defined > > Numpy version: 1.0.1.dev3436 > Scipy version: 0.5.2.dev2319 > Can you help me out or point me to someone else who can? Thanks in > advance. Thank you for your interest in pyem. First, may I request that next time, you set a ticket to scipy tracker ? It is easier to track things, and everybody has access to it. Now, there are actually 3 problems in the code ! A syntax error, as you spot; the example is also wrong because I inverted d and k: GM.__init__ requires first the number of dimension, then the number of components... The error should have been obvious, if there wasn't a third error, that is an error in the raised exception string ! I will update the code + doc accordingly, everything should be available in svn in a few minutes, cheers, David

1 0

Problem converting from scipy_distutils to numpy.distutils
by Berthold Höllmann Nov. 15, 2006

Nov. 15, 2006

Hello, I maintain a larger project where I utilize scipy_distutils for quite some time. Now I find the time to convert the project to numpy from Numeric. I want to be able to use different versions of the INTEL Fortran compiler under Linux easily, so in my path I have only script named like ifort90 or ifort91 that call the correct compiler. So To compile my project I use a command line like: python setup.py config_fc --fcompiler=intel --f90exec=ifort91 --f77exec=ifort91 build It used to work with scipy_distutils, but now I get: ... customize UnixCCompiler using build_ext warning: build_ext: fcompiler=intel is not available. ... Traceback (most recent call last): File "setup.py", line 201, in <module> ext_modules=extInfo.exts) File "/usr/software/gltools/python/Python-2.5/lib/python2.5/site-packages/numpy-1.0-py2.5-linux-i686.egg/numpy/distutils/core.py", line 174, in setup return old_setup(**new_attr) File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/core.py", line 151, in setup dist.run_commands() File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/dist.py", line 974, in run_commands self.run_command(cmd) File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/command/build.py", line 112, in run self.run_command(cmd_name) File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/cmd.py", line 333, in run_command self.distribution.run_command(command) File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/usr/software/gltools/python/Python-2.5/lib/python2.5/site-packages/numpy-1.0-py2.5-linux-i686.egg/numpy/distutils/command/build_ext.py", line 121, in run self.build_extensions() File "/usr/local/gltools/python/Python-2.5/lib/python2.5/distutils/command/build_ext.py", line 407, in build_extensions self.build_extension(ext) File "/usr/software/gltools/python/Python-2.5/lib/python2.5/site-packages/numpy-1.0-py2.5-linux-i686.egg/numpy/distutils/command/build_ext.py", line 312, in build_extension link = self.fcompiler.link_shared_object AttributeError: 'NoneType' object has no attribute 'link_shared_object' It seems I need a command line switch for the linker as well? 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.

1 0

scipy.stats.sem is wrong
by Roman Bertle Nov. 15, 2006

Nov. 15, 2006

Hello, i think scipy.stats.sem is wrong. It gives the same result as scipy.stats.stderr (using N-1 and not N), whereas scipy.stats.tsem uses N and gives the correct result. I have attached a patch correcting this. Related to this, i wonder why there are so many related functions in scipy.stats doing the same, but in a slightly different way. E.g. there are nanstd, std, tstd, some use numpy.std, some not, some take an axis argument, some not. And there is samplestd and samplevar, but sampleerr is called sem instead. Shouldn't these functions be unified somehow? Regards, Roman ------------------------- diff -rud python-scipy-0.5.1/Lib/stats/stats.py python-scipy-0.5.1-new/Lib/stats/stats.py --- python-scipy-0.5.1/Lib/stats/stats.py 2006-08-29 11:58:37.000000000 +0200 +++ python-scipy-0.5.1-new/Lib/stats/stats.py 2006-11-15 12:18:23.000000000 +0100 @@ -1166,9 +1166,7 @@ integer (the axis over which to operate) """ a, axis = _chk_asarray(a, axis) - n = a.shape[axis] - s = samplestd(a,axis) / sqrt(n-1) - return s + return samplestd(a,axis) / float(sqrt(a.shape[axis])) def z(a, score): diff -rud python-scipy-0.5.1/Lib/stats/tests/test_stats.py python-scipy-0.5.1-new/Lib/stats/tests/test_stats.py --- python-scipy-0.5.1/Lib/stats/tests/test_stats.py 2006-08-29 11:58:37.000000000 +0200 +++ python-scipy-0.5.1-new/Lib/stats/tests/test_stats.py 2006-11-15 12:11:29.000000000 +0100 @@ -740,15 +740,16 @@ ## assert_approx_equal(y,0.775177399) y = scipy.stats.stderr(self.testcase) assert_approx_equal(y,0.6454972244) + def check_sem(self): """ this is not in R, so used - sqrt(var(testcase)*3/4)/sqrt(3) + sqrt(samplevar(testcase))/sqrt(4) """ #y = scipy.stats.sem(self.shoes[0]) #assert_approx_equal(y,0.775177399) y = scipy.stats.sem(self.testcase) - assert_approx_equal(y,0.6454972244) + assert_approx_equal(y,0.5590169944) def check_z(self): """ -------------------------

2 1

LAPACK 3.1
by Nils Wagner Nov. 13, 2006

Nov. 13, 2006

FWIW http://www.netlib.org/lapack/lapack-3.1.0.changes Nils

1 0

banded generalized eigenvalue problems
by Jonathan Taylor Nov. 10, 2006

Nov. 10, 2006

Earlier message bounced: wrong email address: ------------------ Hi, I need to solve a banded generalized eigenvalue problem and was going to try to mimic the code in generic_flapack.pyf to generate a wrapper for the appropriate lapack function: dsbgv (and its name variants). Is this the recommended way of using extra lapack functionality that is not presently in scipy? If not, any other suggestions? Thanks, Jonathan

3 2

How to speed up the computation of triple integrals
by Nils Wagner Nov. 7, 2006

Nov. 7, 2006

Hi all, Is there a way to speed up the computation of triple integrals ? Attached is a small script illustrating the task. Any pointer would be appreciated. Nils

2 2

scipy.org serving up Internal Server Errors?
by Albert Strasheim Nov. 6, 2006

Nov. 6, 2006

Hello all scipy.org seems to be serving up Internal Server Error pages intermittently. Alternatively, trying to access something like http://projects.scipy.org/mpi4py/timeline simply hangs after the connection is established. Enthought folks? Cheers, Albert

3 2