Mailman 3 January 2006 - SciPy-User

Matrix Multiply for > 2 dim matricies?
by Gregory Novak Jan. 31, 2006

Jan. 31, 2006

I'm somewhat mystified by the behavior of the matrix multiply routine for > 2d matricies. The back story: at every point on a 3d grid, I have a vector and I would like to apply a matrix transformation to all of them. I had thought that by feeding the matrix and the whole 4d array to matrixmultiply, I'd be able to quickly do this. However, I had trouble figuring out how to correctly interpret the output from matrixmultiply. After exhaustively feeding 2x2x2 arrays to matrix multiply to see if I could tease out the correct behavior, I found that A = matrixmultiply(B,C) means that: A_ij0k = B_ijn C_0nk But when the first index of C is 1, I don't see what's going on. For example, it seems that: A_0010 = B_000 * C_010 + B_001 * C_100 Ie, for either: B = zeros((2,2,2)); C = zeros((2,2,2)) B[0,0,0] = 1 C[0,1,0] = 1 A = matrixmultiply(B,C) or: B = zeros((2,2,2)); C = zeros((2,2,2)) B[0,0,1] = 1 C[1,0,0] = 1 A = matrixmultiply(B,C) A ends up being the same, namely: [[[[0,0,] [1,0,]] [[0,0,] [0,0,]]] [[[0,0,] [0,0,]] [[0,0,] [0,0,]]]] So, I know that matrixmultiply multiplies a bunch of things and adds them together, but I'm having trouble deducing exactly _what_ it does in this case. Thanks, Greg

2 1

numpy and f2py
by Ryan Krauss Jan. 31, 2006

Jan. 31, 2006

Is there a secret to making f2py and the new numpy play together? Here is the begining of a fortran file that worked great with old scipy. I just tried installing numpy/new scipy and it no longer works - the output is a vector of nan's (I have not re-installed f2py since installing new scipy). subroutine bodevect(svect,ucv,outvect,n1,n2) Cf2py integer intent(hide),depend(svect) :: n1 = len(svect) Cf2py integer intent(hide),depend(ucv) :: n2 = len(ucv) Cf2py intent(in) svect, ucv Cf2py intent(out) outvect integer n1, n2, i double complex svect(n1),outvect(n1), bode double precision ucv(n2) DO i=1,n1 outvect(i) = bode(svect(i),ucv,n1) ENDDO END double complex function zcosh(z) double complex z zcosh = 0.5*(exp(z)+exp(-z)) RETURN END double complex function zsinh(z) double complex z zsinh = 0.5*(exp(z)-exp(-z)) RETURN END double complex function bode(s,ucv,n) double complex s double complex zsinh, zcosh Cf2py intent(in) s, ucv Cf2py intent(out) bode Cf2py integer intent(hide),depend(ucv) :: n = len(ucv) integer n double precision ucv(n) double precision kbase, cbase, mubeam, EIbeam, Lbeam, rl0, Ll0, @ ml0, Il0, kj1, cj1, rl1, Ll1, ml1, Il1, Kact, tauact, kj2, @ cj2, rl2, Ll2, ml2, Il2, kj3, cj3, rl36, Ll36, ml36, Il36, @ gainbode1, abeam, gainbode0 double complex c1beam, c2beam, c3beam, c4beam, betabeam, a_1, a_2, @ a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_10, a_11, a_12, a_13, @ a_14, a_15, a_16, a_17, a_18, a_19, a_20, a_21, a_22, a_23, @ a_24, a_25, a_26, a_27, a_28, a_29, a_30, a_31, a_32, a_33, @ a_34, a_35, a_36, a_37, a_38, a_39, a_40, a_41, a_42, a_43, @ a_44, a_45, a_46, a_47, a_48, a_49, a_50, a_51, a_52, a_53, @ a_54, a_55, a_56, a_57, a_58, a_59, a_60, a_61, a_62, a_63, @ a_64, a_65, a_66, a_67, a_68, a_69, a_70, a_71, a_72, a_73, @ a_74, a_75, a_76, a_77, a_78, a_79, a_80, a_81, a_82, a_83, @ a_84, a_85, a_86, a_87, a_88, a_89, a_90, a_91 kbase=ucv(1) Thanks, Ryan

3 5

What happened to Numeric Python EM Project?
by Frederick Ross Jan. 31, 2006

Jan. 31, 2006

On the website there's a link to the Numeric Python EM Project, which looks really interesting from the description...but it appears to have vanished from the face of the earth. www.pythonemproject.com isn't responding, and Rob Lytle's home page doesn't appear to have a link to it. Does anyone know where it went? -- Frederick Ross Graduate Fellow The Rockefeller University

1 0

Centroid Calculation
by Brian Cole Jan. 30, 2006

Jan. 30, 2006

I'm a NumPy newbie. What is the NumPy way of doing this? n=0 centroid_x=0 centroid_y=0 centroid_z=0 for x, y, z in catesian_coords: centroid_x+=x centroid_y+=y centroid_z+=z n+=1 centroid_x/=n centroid_y/=n centroid_z/=n Or am I entirely missing the point of NumPy? Thanks, Brian

2 1

SciPy migration underway
by Joe Cooper Jan. 30, 2006

Jan. 30, 2006

Hi all, As many of you know, there's been a lot of work going on on various aspects of the SciPy website--Andrew Straw and Travis Oliphant and many others have been working hard on the new Moin-based website, and it's looking really great. It has already proven to be a success in all of the ways that the Plone site has always failed (i.e. making it easy for people to get involved and contribute), and everyone concerned would now like it to be the official First Contact for anyone coming to SciPy. To that end tonight I will be making the final changes that moves the website, mailing lists and old archives and user lists, and everything else to the new server and the front page of scipy.org will be the new Moin wiki (http://new.scipy.org/Wiki for those who haven't seen it, but really want to see it even before the move takes place this afternoon). I'm putting forth an effort not to break anything in the process, but that's hoping for a lot when so many services and complicated redirects, proxy rules, and assorted stuff is involved. So, if you find any SciPy.org service is broken between now and tomorrow morning, don't be too surprised. If anything is still broken in the morning, please let me know about it, as it probably means I didn't notice. Anyway, just a heads up about the impending intermittent downtime of some or all of our services on SciPy.org during the next few hours.

2 1

cubic splines in scipy
by Andrew Straw Jan. 30, 2006

Jan. 30, 2006

Hi, (Warning: non math guru alert.) I'm looking to do 1D interpolation using cubic splines: http://mathworld.wolfram.com/CubicSpline.html That page has a 2D picture, but the math is for 1D. The important bit for me is that interpolated curve passes through the control points and the whole thing is continuously differentiable. I can see a few hints toward what I want to do in scipy: scipy.signal.cspline1d and various spline routines in scipy.interpolate. Unfortunately, cspline1d doesn't seem to be tied to anything I can use to interpolate with and the things in scipy.interpolate don't seem to be cubic splines -- I desire my interpolated trajectory to pass through the control points. Can anyone point me to anything in scipy that allows cubic spline interpolation, including the interpolated values passing through the control points? According to the MathWorld page cited above, I thought the attached code would work (not optimized for speed, obviously), but it doesn't interpolate well. I guess from this that scipy.signal.cspline1d isn't doing exactly what I think it's doing. Or perhaps my code is wrong? (For what it's worth, I need to interpolate discretely sampled data from a discrete-time simulation that I then want to use as forcing functions for use in a scipy ODE solver. So, my discretely sampled data contain no noise and I don't want to smooth them, I'm just looking for intermediate values that are continuously differentiable.) Thanks in advance for any help, Andrew

3 5

Re: [SciPy-user] Confusing BLAS/LAPACK situation on Mac OS X
by Brian Granger Jan. 29, 2006

Jan. 29, 2006

On Jan 28, 2006, at 1:06 PM, Robert Kern wrote: Brian Granger wrote: Thus, It looks like the Apple supplied cblas and clapack libraries are never used, even though they are found at build time. What is the current status of using Apple's cblas and clapack? Can I use them from within scipy/numpy? How? AFAICT, Accelerate.framework does not include CLAPACK versions of their LAPACK functions, although they do include CBLAS versions of their BLAS functions. I've never looked into how scipy determines the existince of the CBLAS routines because I've always used ATLAS since Accelerate.framework does not contain the CLAPACK versions. What about: /System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/ Versions/Current/libLAPACK.dylib and Headers/clapack.h Also, Apple's docs on Accelerate claim that they include clapack. Are these incomplete/unusable for some reason? It seems silly to build ATLAS if it is already there. However, it looks like the culprit is line 82 in scipy/Lib/linalg/setup.py. ... 82 if name[0]=='c' and atlas_version is None and newer(__file__,target): f = open(target,'w') f.write('python module '+name+'\n') f.write('usercode void empty_module(void) {}\n') ... I don't think atlas_version gets set when the ATLAS is Accelerate.framework. Maybe I will play around with this, but I am not very familiar with how libalg is setup, so I don't know if I will get anywhere. -- Robert Kern robert.kern(a)gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter _______________________________________________ SciPy-user mailing list SciPy-user(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-user Brian Granger, Ph.D. Assistant Professor of Physics Santa Clara University bgranger(a)scu.edu Phone: 408-551-1891 Fax: 408-554-6965 This message scanned for viruses and SPAM by GWGuardian at SCU (MGW1)

4 5

Re: [SciPy-user] Confusing BLAS/LAPACK situation on Mac OS X
by Brian Granger Jan. 28, 2006

Jan. 28, 2006

Brian Granger, Ph.D. Assistant Professor of Physics Santa Clara University bgranger(a)scu.edu Phone: 408-551-1891 Fax: 408-554-6965 >>> oliphant.travis(a)ieee.org 01/28/06 3:23 PM >>> Arnd Baecker wrote: >With numpy > mat = zeros((Ny, Nx), fortran=1) >should do the job (I have not tested this yet). > >If all the above is correct, then the solution I would use >for myself is to set `fortran=1` for all arrays which >will be used by some LAPACK routine. >Then no unnecessary copies of (presumably large) arrays will take >place on any machine and one could stop worrying about flapack >vs. clapack ;-). > >Not sure if that is the solution for everyone - so I am happy >to learn about any drawbacks.... > > This is the main purpose of the Fortran-order arrays in NumPy --- optimize interfaces to Fortran-written packages. Right now, C-contiguous arrays still have a "special place" because several algorithms require C-contiguous arrays in order to work and will make copies of Fortran-order arrays as needed. There may still be some issues with these Fortran-order arrays especially regarding un-needed copying. With Numeric, f2py did an intelligent job of deciding whether or not to copy. Most of this is unneeded now because the FORTRAN flag on the NumPy array is kept up-to-date with the striding information on the NumPy array so that you just need to look at that flag to determine if the array is in Fortran-order or not. -Travis So it sounds like on Mac OS X, I should use fortran=1 in my numpy arrays if I will be calling things like linalg.eig. That way no copies are made when they are passed to the Accelerate framework's versions of LAPACK which take column major arrays. It is really nice to be able to pick the storage format like this. Brian _______________________________________________ SciPy-user mailing list SciPy-user(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-user This message scanned for viruses and SPAM at SCU (MGW2)

1 0

Re: [SciPy-user] Confusing BLAS/LAPACK situation on Mac OS X
by Brian Granger Jan. 28, 2006

Jan. 28, 2006

Brian Granger, Ph.D. Assistant Professor of Physics Santa Clara University bgranger(a)scu.edu Phone: 408-551-1891 Fax: 408-554-6965 >>> robert.kern(a)gmail.com 01/28/06 1:43 PM >>> Brian Granger wrote: > What about: > > /System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/ > > Versions/Current/libLAPACK.dylib > > and > > Headers/clapack.h > > Also, Apple's docs on Accelerate claim that they include clapack. > > Are these incomplete/unusable for some reason? It seems silly to build ATLAS if it is already there. Look at the symbols in libBLAS.dylib and libLAPACK.dylib. You can see the CBLAS versions in libBLAS.dylib, but there are only the FORTRAN versions in libLAPACK.dylib. You are right that the docs claim (correctly) that Accelerate.framework contains CLAPACK. However, I goofed and referred to things incorrectly. CLAPACK is simply an implementation of FORTRAN LAPACK in C. It uses FORTRAN column-major arrays. The flapack module links against these. However, ATLAS also provides row-major versions of some(?) LAPACK functions to match the row-major CBLAS routines. The clapack module wraps these. These are not included in Accelerate.framework, although the row-major CBLAS versions are. I have been reading the header files in the Accelerate framework to get a better understanding of all this. Apple's implementation of BLAS is the usual CBLAS, which can work with matrices in either row or column major storage. The storage order is determined at runtime by an argument passed to the functions. Are these wrapped and available in scipy? If so, then where? The LAPACK implmentation, as you assert, includes only the fortran style functions that require columnn major arrays. So even though these are called "CLAPACK" by Apple, scipy refers to them as "FLAPACK" because they are not the row major versions that ATLAS exposes. Is this correct? Thanks Brian -- Robert Kern robert.kern(a)gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter _______________________________________________ SciPy-user mailing list SciPy-user(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-user This message scanned for viruses and SPAM at SCU (MGW2)

1 0

Linux Installation Instructions
by John Byrnes Jan. 28, 2006

Jan. 28, 2006

Is anyone else unable to access the Linux install instructions at http://pong.tamu.edu/tiki/tiki-view_blog_post.php?blogId=6&postId=97 I get a 403 - Forbidden Page. Regards, John -- Well done is better than well said. -- Benjamin Franklin

2 1