Mailman 3 November 2009 - SciPy-User

Re: [SciPy-User] Distributed computing: running embarrassingly parallel (python/c++) codes over a cluster
by Luis Pedro Coelho Nov. 12, 2009

Nov. 12, 2009

Rohit Garg wrote: > I have an embarrassingly parallel problem, very nicely suited to > parallelization. I have lots of those :) > My only constraint is that it should be able to run a python extension > (c++) with minimum of fuss. I want to minimize the headaches involved > with setting up/writing the boilerplate code. Which > framework/approach/library would you recommend? My own: It's called jug. See http://luispedro.org/software/jug ( Or download the code from github: http://github.com/luispedro/jug ) * It works with any set of processors that can either share a filesystem (plays well with NFS, but can be slow) or a connection to a redis database (which is very easy to set up and is probably as fast as any other approach if everyone is on the same processor). A major advantage is that you write mostly Python (and not something funny looking). For example, here's what a programme with that framework would look like: @TaskGenerator def preprocess(input): ... @TaskGenerator def compute(input, param): ... @TaskGenerator def collect(inputs): ... results = [] for input in glob('*.in'): intermediate = preprocess(input) results.append(compute(intermediate, param)) final = collect(results) The only step that's different w.r.t. to the linear version is adding the TaskGenerator decorator, which changes a call of preprocess(input) into Task(preprocess, input). Jug handles everything else. I have been using this now for almost year for all my research work and it works very well for me. HTH, Luis

1 0

ANN: scikits.image v0.2
by Stéfan van der Walt Nov. 12, 2009

Nov. 12, 2009

I'm glad to announce the second release of `scikits.image`, a collection of image processing routines for SciPy. On top of bug-fixes and improved documentation, the following changes/additions were made: - A new IO plugin infrastructure so that commands like 'imshow' are available via multiple backends (PIL, matplotlib, QT4, etc.) - ImageCollections (for cached loading of multiple images) and MultiImage (for working with multi-layered images) - More complete OpenCV wrappers - A graphical image viewer (also installed as a script `scivi`), that allows colour adjustments - Shortest path algorithm For version 0.3, we aim to - Incorporate some of the code offered by the Broad institute - Implement acquisition (grabbing images from cameras) and intrinsic camera calibration - Add real time video and camera display with processing - Improve filtering code - Add morphological operations More information is available at: http://stefanv.github.com/scikits.image/ Regards Stéfan

1 0

vectorplot scikit code patch
by George Nurser Nov. 12, 2009

Nov. 12, 2009

Hi, Not sure where to post this. The lic_demo.py and lic_efield_demo.py scripts in the vectorplot scikit fail for me with Traceback (most recent call last): File "lic_efield_demo.py", line 55, in <module> plt.figimage(image) File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/pyplot.py", line 404, in figimage sci(ret) File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/pyplot.py", line 160, in sci gca()._sci(im) File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/matplotlib/axes.py", line 1338, in _sci "Argument must be an image, collection, or ContourSet in this Axes") ValueError: Argument must be an image, collection, or ContourSet in this Axes I got both scripts to work by in each of them replacing plt.clf() plt.axis('off') plt.figimage(image) by fig = plt.figure() plt.clf() plt.axis('off') fig.figimage(image) --George.

1 0

loading mat file in scipy
by Alireza Alemi-Neissi Nov. 12, 2009

Nov. 12, 2009

Dear All, I do not know if it is the right place to ask this question. If not, please let me know where I can post it. This the problem: I have recently tried to use scipy.io.loadmat() to load a mat file to python. The mat file is a complicated nested structure in matlab. It turned out that it takes a long time (around 5 min) to load this data the size of which is around 120MB. On the other hand it take only a couple of seconds to load it in Matlab. What is wrong with loadmat in scipy? Thanks, Alireza ---------------------------------------------------------------- SISSA Webmail https://webmail.sissa.it/ Powered by SquirrelMail http://www.squirrelmail.org/

5 9

Use of MPI in extension modules
by Stefan Seefeld Nov. 11, 2009

Nov. 11, 2009

Hello, I have a rather basic question about using (C++) extension modules with ipython. Sorry if this is the wrong list for this. I'm working on a signal & image processing library that uses MPI internally. I'd like to provide a Python interface to it, so I can integrate it into SciPy. With 'normal' Python this all works nicely. Just recently I have started to consider parallelism, i.e. I want to use the library's internal parallelism, by running it with ipython in parallel. My assumption was that all the engines started via 'ipcluster mpiexec ..." would already have MPI_Init called, and thus, my extension modules would merely share the global MPI state with the Python interpreter. That doesn't seem to be the case, as I either see all my module instances report rank 0, or, if I don't call MPI_Init, get a failure on the first MPI call I do. Can anybody help ? Do I need to initialize MPI myself in my extension module ? Any pointers are highly appreciated. Thanks, Stefan -- ...ich hab' noch einen Koffer in Berlin...

4 5

the skellam distribution
by Ernest Adrogué Nov. 10, 2009

Nov. 10, 2009

Hi, In case somebody is interested, or you want to include it in scipy. I used these specs here from the R package: cran.r-project.org/web/packages/skellam/skellam.pdf Note that I am no statician, somebody who knows what he's doing (as opposed to me ;) should verify it's correct. import numpy import scipy.stats.distributions # Skellam distribution ncx2 = scipy.stats.distributions.ncx2 class skellam_gen(scipy.stats.distributions.rv_discrete): def _pmf(self, x, mu1, mu2): if x < 0: px = ncx2.pdf(2*mu2, 2*(1-x), 2*mu1)*2 else: px = ncx2.pdf(2*mu1, 2*(x+1), 2*mu2)*2 return px def _cdf(self, x, mu1, mu2): x = numpy.floor(x) if x < 0: px = ncx2.cdf(2*mu2, x*(-2), 2*mu1) else: px = 1-ncx2.cdf(2*mu1, 2*(x+1), 2*mu2) return px def _stats(self, mu1, mu2): mean = mu1 - mu2 var = mu1 + mu2 g1 = (mu1 - mu2) / numpy.sqrt((mu1 + mu2)**3) g2 = 1 / (mu1 + mu2) return mean, var, g1, g2 skellam = skellam_gen(a=-numpy.inf, name="skellam", longname='A Skellam', shapes="mu1,mu2", extradoc="") Bye. -- Ernest

4 13

Saving record arrays to tab formatted txt
by David Baddeley Nov. 10, 2009

Nov. 10, 2009

Hi all, does anyone know of an easy way to save record arrays as tab formatted txt? numpy.savetxt doesn't do the trick. I've got a record array with nested records, and mixed data types - dtype is as follows: dtype([('tIndex', '<i4'), ('fitResults', [('A', '<f4'), ('x0', '<f4'), ('y0', '<f4'), ('sigma', '<f4'), ('background', '<f4'), ('bx', '<f4'), ('by', '<f4')]), ('fitError', [('A', '<f4'), ('x0', '<f4'), ('y0', '<f4'), ('sigma', '<f4'), ('background', '<f4'), ('bx', '<f4'), ('by', '<f4')]), ('resultCode', '<i4'), ('slicesUsed', [('x', [('start', '<i4'), ('stop', '<i4'), ('step', '<i4')]), ('y', [('start', '<i4'), ('stop', '<i4'), ('step', '<i4')]), ('z', [('start', '<i4'), ('stop', '<i4'), ('step', '<i4')])])]) and want to flatten it so that each entry in the table becomes a row in the .txt file. ie: tIndex fitResults_A fitResults_x0 ....\n tIndex fitResults_A fitResults_x0 ....\n etc... If necessary I'll write my own code to format up a string for each line, but thought I'd ask first in case anyone knew of a pre-existing solution. It'd ideally be generic as I'm also likely to want to use variations on the data type. many thanks, David

2 1

Plotting across the Prime Meridian
by Bruce Ford Nov. 10, 2009

Nov. 10, 2009

In the code below (which is an extraction from a larger set of file), I can plot fine across the dateline and elsewhere, however any plots across the PM error out. I'm using OpenLayers to allow users to designate their plot area via a web form by dragging a map area. OpenLayers provide longitudes from 180/-180. The data I'm extracting has longitudes from 0-360 so some manipulation is necessary. I suspect the problem that is causing the error is the line that reads: "grid = grid[slat:nlat+1,wlong:elong+1]" because my wlong is greater than my elong. Does anyone have a suggestion for handling this at the PM. Any suggestions would be appreciated. The error is in the comments of the below script: import matplotlib import matplotlib.pyplot as pyplot #used to build contour and wind barbs plots import matplotlib.colors as pycolors #used to build color schemes for plots import numpy.ma as M #matrix manipulation functions from mpl_toolkits.basemap import Basemap import numpy as np #used to perform simple math functions on data from numpy import * ############################################################## def read_netcdf(filename, hour, param, wlong, elong, nlat, slat,day_time_level): from netCDF4 import Dataset #interprets NetCDF files nc_file = Dataset(filename, mode="r") nlat = nlat+90 slat = slat+90 swh = nc_file.variables['sig_wav_ht'][:] swh = np.squeeze(swh) grid = swh[int(day_time_level)-1, :,: ] grid = M.array(grid) grid = M.masked_where(grid < 0.001, grid) grid = grid[slat:nlat+1,wlong:elong+1] return grid ############################################################### map_res = "i" nlat = 7 slat = -7 #******************************************************* #when the below values for wlong and elong are as such # (across the Prime Meridian) I get the following error: #Traceback (most recent call last): # File "C:\Program Files\Wing\src\debug\tserver\_sandbox.py", line 58, in <module> # File "C:\Python25\Lib\site-packages\numpy\ma\core.py", line 4262, in min # result = self.filled(fill_value).min(axis=axis, out=out).view(type(self)) #ValueError: zero-size array to ufunc.reduce without identity wlong = -12 elong = 8 #******************************************************* # When the below two values are as such, I get the figure I intend #wlong = -89 #elong = -75 day_raw = 1 param = "sig" year1 = 1993 month = "01" vint = 15 para_spacing = 10 merid_spacing = 10 hour = "%02d" % day_raw if wlong < 0: wlong = wlong+360 if elong < 0: elong = elong+360 vtype = "auto" filename = "x:/ww3/NetCDF/daily_means/ww3dm."+str(year1)+str(month)+ ".nc" grid = read_netcdf(filename, hour=hour, param=param, wlong=wlong, elong=elong, nlat=nlat, slat=slat, day_time_level=day_raw) m = Basemap(projection='cyl',resolution=map_res,llcrnrlon=wlong,llcrnrlat=slat,urcrnrlon=elong,urcrnrlat=nlat) x,y = m(*np.meshgrid(range(wlong,elong+1),range(slat,nlat+1))) print wlong print elong if vtype == 'auto': vmin = grid.min() print "<h3>Vmin: ",vmin,"</h3>" vmax = grid.max() print "<h3>Vmax: ",vmax,"</h3>" #vmin = 0 #vmax = 30 elif vtype== "manual": grid = M.masked_outside(grid,vmin,vmax) pyplot.jet() plot = m.contour(x,y,grid,int(vint)-1,linewidths=0.5,colors='k') plot = m.contourf(x,y,grid,int(vint)-1,cmap=pyplot.cm.jet) m.drawcoastlines() #draw coastlines m.drawmapboundary() #draw a line around the map region m.fillcontinents(color='0.8', lake_color=None, ax=None, zorder=None) #fill in continents with color (gray) m.drawparallels(np.arange(-90,90,para_spacing),labels=[1,0,0,0]) #draw parallels m.drawmeridians(np.arange(-180,180,merid_spacing),labels=[0,0,0,1]) #draw meridians pyplot.show() Bruce --------------------------------------- Bruce W. Ford Clear Science, Inc. bruce(a)clearscienceinc.com bruce.w.ford.ctr(a)navy.smil.mil http://www.ClearScienceInc.com Phone/Fax: 904-379-9704 8241 Parkridge Circle N. Jacksonville, FL 32211 Skype: bruce.w.ford Google Talk: fordbw(a)gmail.com

1 0

Contribution to Performance Python
by Ramon Crehuet Nov. 10, 2009

Nov. 10, 2009

Robin, I know it is a wiki, but you need to get permission to modify it. I've tried registering at: http://docs.scipy.org/numpy/accounts/login but I always get an Authentication failed error. However if I try to register again, it complains that my username is already in use (it wasn't the first time I registered, so that is me! :-) ) My username is rcrehuet, and the address rcsqtc_at_iqac.csic.es Cheers, Ramon PS. David, I am willing to contribute to the cookbook. On Fri, Nov 6, 2009 at 10:53 AM, Ramon Crehuet <rcsqtc(a)iqac.csic.es> wrote: > > If this is interesting to the community, who should I contact to have > > this included in the scipy web page? Hi, It's a wiki, so I think you should be able to register an account and modify the page yourself. I'd certainly support the addition of the new fortran versions (although perhaps one is sufficient since they both seem to perform very closely) and perhaps updated timings section. Cheers Robin

2 1

Contribution to Performance Python
by Ramon Crehuet Nov. 10, 2009

Nov. 10, 2009

Hi all, After reading the "Performance Python" page at: http://www.scipy.org/PerformancePython?action=show I thought some code with Fortran 90/95 was missing, in partuclar considering its useful array features. So I have written a couple of examples with Fortran 90 arrays and with the Fortran 95 forall construct. Both are nicer (to me :-) ) than the FORTRAN77 loops and also faster! In my PC a 1000x1000 array gives: Doing 100 iterations on a 1000x1000 grid numeric took 5.86 seconds fortran77 took 3.53 seconds fortran90-arrays took 1.58 seconds fortran95-forall took 1.58 seconds slow (1 iteration) took 9.13 seconds 100 iterations should take about 913.000000 seconds If this is interesting to the community, who should I contact to have this included in the scipy web page? Cheers, Ramon This are the two new subroutines. I can send the modified laplace.py to whoever wants it. ****************************************************** ! File flaplace90_arrays.f90 subroutine timestep(u,n,m,dx,dy,error) implicit none real (kind=8), dimension(0:n-1,0:m-1), intent(inout):: u real (kind=8), intent(in) :: dx,dy real (kind=8), intent(out) :: error integer, intent(in) :: n,m real (kind=8), dimension(0:n-1,0:m-1) :: diff real (kind=8) :: dx2,dy2,dnr_inv !f2py intent(in) :: dx,dy !f2py intent(in,out) :: u !f2py intent(out) :: error !f2py intent(hide) :: n,m dx2 = dx*dx dy2 = dy*dy dnr_inv = 0.5d0 / (dx2+dy2) diff=u u(1:n-2, 1:m-2) = ((u(0:n-3, 1:m-2) + u(2:n-1, 1:m-2))*dy2 + & (u(1:n-2,0:m-3) + u(1:n-2, 2:m-1))*dx2)*dnr_inv error=sqrt(sum((u-diff)**2)) end subroutine ****************************************************** ! File flaplace95_forall.f90 subroutine timestep(u,n,m,dx,dy,error) implicit none real (kind=8), dimension(0:n-1,0:m-1), intent(inout):: u real (kind=8), intent(in) :: dx,dy real (kind=8), intent(out) :: error integer, intent(in) :: n,m real (kind=8), dimension(0:n-1,0:m-1) :: diff real (kind=8) :: dx2,dy2,dnr_inv integer :: i,j !f2py intent(in) :: dx,dy !f2py intent(in,out) :: u !f2py intent(out) :: error !f2py intent(hide) :: n,m dx2 = dx*dx dy2 = dy*dy dnr_inv = 0.5d0 / (dx2+dy2) diff=u forall (i=1:n-2,j=1:m-2) u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2+(u(i,j-1) + u(i,j+1))*dx2)*dnr_inv end forall error=sqrt(sum((u-diff)**2)) end subroutine ******************************************************

3 2