Mailman 3 March 2010 - SciPy-User

Lagrange Multipliers in optimize.slsqp
by Guilherme P. de Freitas 10 Mar '10

10 Mar '10

Hi all, Would it be possible for optimize.slsqp to return the Lagrange multipliers of the constraints? I don't know Fortran, but it seems that they are available in the original Fortran code: http://projects.scipy.org/scipy/browser/trunk/scipy/optimize/slsqp/slsqp_op… Thanks, Guilherme -- Guilherme P. de Freitas http://www.gpfreitas.com

4 7

<transferred from Scipy-Dev>Re: Seeking help/ advice for applying functions
by eat 10 Mar '10

10 Mar '10

> > On 9 March 2010 09:59, eat <e.antero.tammi <at> gmail.com> wrote: > > Robert Kern <robert.kern <at> gmail.com> writes: > > > >> > >> Your example is not very clear. Can you write a less cryptic one with > >> informative variable names and perhaps some comments about what each > >> part is doing? > >> > > Anne Archibald <peridot.faceted <at> gmail.com> writes: Hi, I moved this thread from scipy-dev to scipy-user. First of all, Thanks Anne you clarified a lot. > > """ > > Hi, > > > > I have tried to clarify my code. First part is the relevant one, the rest is > > there just to provide some context to run some tests. > > The short answer is that no, there's no way to optimize what you're doing. > > The long answer is: when numpy and scipy are fast, they are fast > because they avoid running python code: if you add two arrays, there's > only one line of python code, and all the work is done by loops > written in C. If your code is calling many different python functions, > well, since they're python functions, to apply them at all you must > necessarily execute python code. There goes any potential speed > advantage. (There may be a convenience advantage; if so, you can look > into using np.vectorize, which is just a wrapper around a python loop, > but is convenient.) Not only that I'm new to Scipy/ Numpy but I'm new to python aswell. I think I assumed that once python functions get compiled they'll would be (almost) as efficient than builtins. > > That said, I assume you are considering numpy/scipy because you have > arrays of thousands or more. It also seems unlikely that you actually > have thousands of different functions (that's an awful lot of source > code!). So if your "different" functions are actually just a handful > (or fewer) pieces of actual code, and you are getting your thousands > of functions by wrapping them up with parameters and local variables, > well, now there are possibilities. Exactly what possibilities depend > on what your functions look like - which is one reason Robert Kern > asked you to clarify your code - but they all boil down to rearranging > the problem so that it goes back to "few functions, much data", then > writing the functions in such a way that you can use numpy to apply > them to thousands or millions of data points at once. Yes, indeed I don't have so many different 'base' functions, but just configured (parametrized) different ways. However there are situations when it's advantage to be able to treat your fuctions as 'black boxes'. For example when parametrization is based on optimization. Yes, I'll also recognize that my original question was too generic. I'll try to be more specific next time. > > > Also originally I should have posted this to Scipy-User list. Would it be > > more appropriate to continue the discussion there? > > Probably. > > Anne Thanks a lot, eat > > > Regards, > > eat > > """ > > > > import numpy as np > > > > ## relevant part are in scipy-dev

4 12

return "full_output" or how to stop throwing away already calculated results
by josef.pktd＠gmail.com 10 Mar '10

10 Mar '10

scipy functions throw away a lot of intermediate results examples fmin_slsqp ticket:1129 scipy.linalg.lstsq throws away residuals : if n<m: x1 = x[:n] if rank==n: resids = sum(x[n:]**2,axis=0) x = x1 return x,resids,rank,s in statsmodels, I am switching more to the pattern when full_output=True (or similar keyword) for a function, then all intermediate results are returned, attached to a generic class instance ( or could be a struct in matlab or a bunch for Gael,...) Is it possible to start the same pattern in scipy? E.g. for optimizers like slsqp, if we want to change the returns then each time it breaks the API or we have to increase full_output= 1,2,3 ( like nargout in matlab) (ugly) I saw linalg.lstsq by chance this morning, but I assume there is additional useful information hidden in other functions. Any opinion ? Josef

3 6

fromimage/imread segfaults on my images
by Daniel Lenski 09 Mar '10

09 Mar '10

Hi all, I'm using SciPy to process some sparse binary images with a resolution of about 5000x5000 pixels. I try to load images like this: im = Image.open('/tmp/foo.pbm') print "Image loaded, size is %s and mode is %s." % (im.size, im.mode) arr = fromimage(im) or arr = imread('/tmp/foo.pbm') In either case, I get a segfault. The segfault is definitely in the SciPy code, rather than the PIL code, since PIL's Image.open runs fine. I ran strace, and there's an attempt to mmap() a region as large as the image, one byte per pixel, immediately before the segfault: mmap(NULL, 21463040, PROT_READ|PROT_WRITE, MAP_PRIVATE|MAP_ANONYMOUS, -1, 0) = 0x7f278b5be000 The funny thing is... the imread() implementation in pylab *does* work, but it converts all images to RGB mode and reverses them vertically for no apparent reason. So I'd rather not use it. Is this a known bug in SciPy? Any workarounds/fixes? I don't know any other good, reliable way to get a big image into an array :-( Dan

2 5

First attempt to import scipy.stats fails
by PHobson＠Geosyntec.com 09 Mar '10

09 Mar '10

But the second one works fine -- so it's not really an issue that bothers me, but I thought it might interest some folks. Here's an IPython session that demonstrates the issue: Python 2.6.2 (r262:71605, Apr 14 2009, 22:40:02) [MSC v.1500 32 bit (Intel)] ... IPython 0.10 -- An enhanced Interactive Python. ... In [1]: import numpy as np In [2]: import scipy as sp In [3]: import scipy.stats as st --------------------------------------------------------------------------- ValueError Traceback (most recent call last) c:\stuff\utils\py\<ipython console> in <module>() C:\Python26\lib\site-packages\scipy\stats\__init__.py in <module>() 5 from info import __doc__ 6 ----> 7 from stats import * 8 from distributions import * 9 from rv import * C:\Python26\lib\site-packages\scipy\stats\stats.py in <module>() 201 202 #import scipy.stats #is this a circular import ? --> 203 from morestats import find_repeats #is only reference to scipy.stats 204 import distributions 205 C:\Python26\lib\site-packages\scipy\stats\morestats.py in <module>() 5 import statlib 6 import stats ----> 7 import distributions 8 from numpy import isscalar, r_, log, sum, around, unique, asarray 9 from numpy import zeros, arange, sort, amin, amax, any, where, \ C:\Python26\lib\site-packages\scipy\stats\distributions.py in <module>() 25 from scipy.special import gammaln as gamln 26 from copy import copy ---> 27 import vonmises_cython 28 import textwrap 29 c:\stuff\utils\py\numpy.pxd in scipy.stats.vonmises_cython (scipy\stats\vonmises_cython.c:2939)() ValueError: numpy.dtype does not appear to be the correct type object In [4]: import scipy.stats as st In [5]: np.__version__ Out[5]: '1.4.0' In [6]: sp.__version__ Out[6]: '0.7.1' Paul M. Hobson Senior Staff Engineer -- Geosyntec Consultants 55 SW Yamhill St, Ste 200 Portland, OR 97204 Phone: 503.222.9518 www.geosyntec.com

2 2

Shared memory using multiprocessing.sharedctypes
by David Baddeley 09 Mar '10

09 Mar '10

I'd like shared memory numpy arrays for some code I'm trying to parallelise. The solution should ideally be cross platform. I remember some earlier discussion about this and, as far as I recall Sturla had a working implementation which was (unfortunately) linux only (please correct me if I'm wrong). Looking a bit further, multiprocessing.sharedctypes seems to have all the stuff needed to handle the low level bits on at least linux and windows, and it seems as though it should be relatively easy to cobble something together on top of this which let you use shared memory numpy arrays relatively transparently. This would have the advantage that all the platform specific nasty stuff is being handled by the multiprocessing module, which as it's now part of the standard library, should hopefully be well maintained & has a good chance of being ported to additional platforms Before I jump in and try writing such a wrapper - essentially a helper function for creating shared numpy arrays (easy) and some pickle handlers to make sure what comes out the other end looks like an array (would need to think/read a bit more, but should be doable) I wanted to make sure that there's not something out there already & that I'm not going to be reinventing the wheel. cheers, David

1 0

Re: [SciPy-User] Compiling scipy.sparse.sparsetools extension throws g++ internal error
by Diederik van Liere 09 Mar '10

09 Mar '10

> Diederik van Liere wrote: >>* *>>* ImportError: /usr/local/lib/python2.6/site-packages/scipy/linalg/clapack.so:*>>* undefined symbol: clapack_sgesv* > what is the output of python setup.py config inside scipy build ? I am not quite sure which output you mean, do you mean the console output or a particular log file? I didn't see any log file in particular. The console output does not give any errors as far as I can tell. Import scipy in python does work. > Also, the output of ldd on clapack.so may be useful. The output of ldd looks suspect (I compared it to other posts): ldd /usr/local/lib/python2.6/site-packages/scipy/lib/lapack/clapack.so libgfortran.so.1 => /usr/lib/libgfortran.so.1 (0x00b96000) libm.so.6 => /lib/libm.so.6 (0x00f41000) libgcc_s.so.1 => /lib/libgcc_s.so.1 (0x00552000) libc.so.6 => /lib/libc.so.6 (0x00110000) /lib/ld-linux.so.2 (0x00e8f000) However, Scipy does find the lapack, blas etc as this snippet from the log shows: atlas_threads_info: Setting PTATLAS=ATLAS libraries lapack_atlas not found in /usr/local/include numpy.distutils.system_info.atlas_threads_info Setting PTATLAS=ATLAS Setting PTATLAS=ATLAS FOUND: libraries = ['lapack', 'lapack', 'f77blas', 'cblas', 'atlas'] library_dirs = ['/usr/local/include'] language = f77 How can I prevent starting a new thread every time? Best, Diederik

3 2

Re: [SciPy-User] scipy+lapack array copy
by Bruce Southey 09 Mar '10

09 Mar '10

On 8. mars 2010, at 10.10, Skipper Seabold wrote: > On Mon, Mar 8, 2010 at 12:47 PM, Paul Anton Letnes > <paul.anton.letnes <at> gmail.com> wrote: >> Hi everyone. >> >> I have a numerical project where I link a fortran module to python using f2py. The fortran module basically sets up the left hand side coefficient matrix of a linear equation system. The right hand side is small, so I set that up in python. So far, everything works great. >> >> Now, when I call scipy-lapack, I am using: >> from scipy.linalg.lapack import flapack as lapack >> lu, piv, x, info = lapack.zgesv(A, b) >> A is created in the fortran part of the program. b is created from numpy.zeros(). Both report as F_CONTIGUOUS by printing A.flags and b.flags. >> >> The problem is that I'd like A to be large (>1GB), but the lapack wrapper (or lapack itself?) copies the array, doubling the memory use. Aside from being unnecessary, this also reduces the practical problem size I can solve by a factor 2, which of course is very annoying. >> >> How can I get around this problem? >> If memory usage is a concern then you can directly form the X.T*X matrix rather than doing the transpose and multiplication. The advantage is that X.T*X has a fixed dimensions regardless of the number of observations. Further, you can exploit symmetry and use half-storage (there are lapack routines). The disadvantage is that is more complex and probably not as fast as doing the multiplication. If numerical stability is a concern then you probably need to use a better solving than LU with partial pivoting anyhow. Bruce

1 0

Re: [SciPy-User] Compiling scipy.sparse.sparsetools extension throws g++ internal error
by Diederik van Liere 08 Mar '10

08 Mar '10

> Diederik van Liere wrote: >>* I need scipy.linalg although I prefer a complete working version of*>>* Scipy because it might the case the sparse module is required by one*>>* of other libraries I am suing.* > Indeed - not building scipy.sparse is clearly a temporary workaround, I > would not advise to do this as a long term solution. Another obvious > thing to try is building the scipy's trunk, or re-generating the c++ > sources from swig. The latter, by generating a different file, may be > able to exercise another codepath in g++, and avoid the crash depending > on where the crash happens in g++. Thanks David, I disable the sparse package and scipy builds fine (I am already using the latest SVN version of Scipy). However,when I import scipy.linalg I receive the following error in Python: ImportError: /usr/local/lib/python2.6/site-packages/scipy/linalg/clapack.so: undefined symbol: clapack_sgesv Which looks very similar to this bug report: http://projects.scipy.org/scipy/ticket/1031 I would like to give you more detailed information but I am not sure what you need / or how to provide that. > Ideally, if you could find-out what part of the code causes the crash > and report the bug, we may be able to implement a workaround in scipy. > cheers, > David PS: I hope this reply does not start a new thread, don't know how to prevent that.

2 1

Compiling scipy.sparse.sparsetools extension throws g++ internal error
by Diederik van Liere 08 Mar '10

08 Mar '10

> Diederik van Liere wrote: >>* error: Command "g++ -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O3*>>* -Wall -fPIC -I/usr/local/lib/python2.6/site-packages/numpy/core/include*>>* -I/usr/local/include/python2.6 -c scipy/sparse/sparsetools/csr_wrap.cxx -o*>>* build/temp.linux-i686-2.6/scipy/sparse/sparsetools/csr_wrap.o" failed with*>>* exit status 1* > What's the version of g++ ? I am using GCC 4.1.2 (32bit) >>* and I really don't know how to fix this. I have tried both Numpy 1.3.0 and*>>* the latest Numpy SVN version but to no avail. I hope someone can point me in*>>* the right direction.* > Changing numpy version is unlikely to fix anything here, as the error is > a crash from the C++ compiler. Okay, good to know. > One solution would be to update gcc (by compiling your own version), or > more simply, if you don't need scipy.sparse to disable it in > scipy/setup.py by commenting the sparse line. I don't think I will need this, this is for sparse matrix manipulation, right? I need scipy.linalg although I prefer a complete working version of Scipy because it might the case the sparse module is required by one of other libraries I am suing. But thanks David for giving me some workarounds. > David

2 1