Mailman 3 August 2009 - SciPy-User

Open Research and Python's role
by Gökhan Sever 13 Aug '09

13 Aug '09

The SciPy09 conference will provide such a great opportunity to bring users and scientists from varying disciplines. I am curious to know whether there is a planned session for a discussion of philosophical outcomes of using an open-source programming language and various other open scientific tools for ones scientific data analysis, visualization, and publishing needs rather than the practical aspects of the subject. Although not a direct Python specific discussion issue, topics such as: - How computational sciences are evolving by the growing open-source communities? - The role of open languages and tools for the science education - Reproducibility, open information, data etc.. (not to mention Popper) You might consider seeing the following resources to get more in-depth information on some of the aforementioned subjects: http://www.openscience.org/blog/?p=269 (The OpenScience Project » What, exactly, is Open Science?) http://hyfen.net/out/writing/2009-07/michael-nielsen-at-science20/ (Doing it in the open: Michael Nielsen at Science 2.0) http://www.linux.com/news/software/applications/32332-open-source-open-rese… http://www.stanford.edu/~vcs/ (Victoria Stodden's Webpage - various talks and papers related to the subject) I still think that these could be an interesting lunch or dinner time talk if not a panel discussion. Gökhan

1 0

vectorizing a function with non scalar arguments
by Luke 13 Aug '09

13 Aug '09

Thanks for the replies! >On 10-Aug-09, at 5:28 PM, Luke wrote: > >> I have a function of the following form: >> def f(q, qd, parameter_list): >> # Unpacking the parameters >> m, g, I11, I22, I33 = parameter_list >> # Unpacking the qdots >> q1, q2, q3, q4, q5, q6 = q >> q1p, q2p, q3p, q4p, q5p, q6p = qd >> ....Some calculations to determine u1,u2,u3,u4,u5,u6.... >> return [u1, u2, u3, u4, u5, u6] > >We'd really need to see "...Some calculations..." but if it's >something that's expressible as simple arithmetic/calls to numpy >ufuncs, then it's probable that you can pull this off. The calculations would involve the basic arithmetic operators (+-*/**) and some trigonmetric functions (sin,cos,tan), and maybe abs, so yes, they could be expressed exclusively as Numpy ufuncs. > >> I would like to vectorize this function so that I could pass it: >> 1) q and qd as a 2-d numpy arrays of shape (n, 6) >> 2) parameter_list as a 1-d numpy array of shape (5,) > >hsplit() and hstack() to unpack/pack the 2D arrays, respectively >(though there are more memory efficient/cache friendly ways, if you >care about speed). I was hoping to code it in such a way that if I could still call the function using 1-d lists, thus making it versatile enough to handle both a single set of q and qd or a set of n q and qd. > >> and it would return a 2-d numpy array of shape (n, 6) >> >> parameter_list is a bunch of constants that don't change, while q and >> qd are the things that are different between say q[n, :] and q[n+1, >> :]. > >> I realize that I can use a for loop to loop from 0 to n-1. >> >> My question is whether there a way to use vectorize (or some other >> decorator), to obtain the behavior? Or is there another nice method >> that people might be able to recommend? > >Again, if you split into 1d arrays and are just doing arithmetic on >them (or you use numpy functions that can take scalar or vector >arguments) then it should "vectorize itself". But it might be helpful >if you have any helper functions of one scalar variable that you need >to evaluate across an entire array of q1's, etc. > >> The documentation for vectorize is very sparse and seems like it is >> geared only towards function which have scalar arguments. > >Yes, 'vectorize' won't do it for you in this case. > Ok. Thanks! >David

1 0

scipy odeint question
by Luke 13 Aug '09

13 Aug '09

I am generating some ODE's that I integrate with scipy.odeint, and everything is working well, but there is something that I would like to do that I'm not sure if odeint can currently handle. The function representing the right hand side of the ODE's I generate are of the standard form that scipy needs, namely: def f(t, x, params): .... return dxdt In the ... there are many (sometimes thousands or tens of thousands of lines) intermediate variables that are used for the final calculation of dxdt. The calculations generally involve the basic arithmetic operators (+-*/**) and some trig functions (sin, cos, tan). There are several other functions of t, x, and params that I would like to evaluate *only at the integration step times specified by the user*, and I would like to be able make use of the fact that the computation of the intermediate variables has already been done for the time step that is to be returned in the x array. Here is an example of what I am trying to acheive: def f(t, x, params): z1 = ... z2 = ... (z1, z2, z3 would be in reality be more like z1, ..., z1000, all functions of t, x, params) ... zn = ... dxdt = [z1+z2**2, z2*z3-z1**3, ...] # The following only to be evaluated at time points in the user specified vector, e.g., t = arange(t_i, t_f, t_s) eqn_1 = f_1(z1,...., zn) eqn_2 = f_2(z1, ..., zn) .... eqn_n = f_n(z1, ..., zn) return dxdt >From a practical point of view, I can easily integrate the equations, determine the state trajectory, then call a function that evaluates the other quantities eqn1, eqn2... *after* numerical integration. But for certain applications, such as realtime control applications, it is useful to avoid repeated calculations (as would happen if I just called a function that took that state trajectory and calculated the eqn1, eqn2). I can see a couple of approaches that might work: 1) Make the time vector a global variable, and use an if statement inside of f(t, x, params) like: if t in t_span: eqn_1 = ... eqn_2 = ... .... eqn_n = ... accumulate_eqns(eqn1, ..., eqn_n) # accumulate_eqns could be either a function passed to f through pararms, or a global function 2) Make eqn_1, ..., eqn_n globals, and only integrate one time step at a time so that they can be stored at the end of each integration step. In this way, odeint would be called in a for loop that would use a time vector that had only two elements (their difference would be the step time specified by the user) I'm not 100% sure the above methods would work because I don't know if globals in the namespace which calls odeint would be viewable to globals declared in the function f. What would be really awesome would be if odeint could hold on to any variables in the f namespace when it reaches a time step that is requested, and then call another function to evaluate the output equations, but that function would have access to the state of all the variables in the f namespace *at the time point* where the state x is to be evaluated. Ideally, the extra equations wouldn't be computed except at the time points which are specified by the time vector specified by the user. Examples of use would be computing the linear/angular momentum of a system of rigid bodies and particles, center of mass of the bodies, total energy (kinetic and potential), and the value of constraints / conserved quantities during the numerical integration of the equations of motion. The equations of motion can be integrated without computing these quantities, but they are very useful for checking that the numerical integration is well behaved. For complicated systems such as spacecraft, these equations can be extremely long, and it can be very computationally advantageous to make use of these intermediate quantities rather than having to be recalculated after the numerical integration is complete. Any ideas? Sincerely, ~Luke

1 0

Re: [SciPy-User] SciPy-User Digest, Vol 72, Issue 32
by Tom Kuiper 13 Aug '09

13 Aug '09

I'm interested in observing (in the United Nations sense of the word) this session. Thursday after the dinner is a bit of a strain since I need to drive back to West LA (30 min) and be up again at 6:30 the next morning. If it's not many people, we could get an early lunch (about 11:30 I think) from the catering truck and find a table in the courtyard of Watson, adjacent to Steele. Tom > Message: 11 > Date: Wed, 12 Aug 2009 19:21:52 -0400 > From: Christopher Hanley <chanley(a)stsci.edu> > Subject: Re: [SciPy-User] Astronomy BoF meeting at SciPy? > To: SciPy Users List <scipy-user(a)scipy.org> > Cc: Christopher Hanley <chanley(a)stsci.edu> > Message-ID: <A95D44F7-1F14-4390-B382-50C618B07894(a)stsci.edu> > Content-Type: text/plain; charset=US-ASCII; format=flowed > > I would prefer to meet Thursday after dinner. > > Chris > > > -- > Christopher Hanley > Senior Systems Software Engineer > Space Telescope Science Institute > 3700 San Martin Drive > Baltimore MD, 21218 > (410) 338-4338 > > On Aug 12, 2009, at 6:14 PM, James Turner wrote: > > >> Hi everyone, >> >> I haven't heard anything regarding an astronomy BoF session since >> Joe's AstroPy post on 26 June regarding astronomy libraries. However, >> at least Joe, Perry and I are interested in talking and I suspect >> there are others (AplPy people?). Up for discussion will be issues >> like co-ordinating library development, as different groups help ramp >> up the functionality available to match IRAF or IDL. >> >> Who else is interested? Do you have existing plans to meet at SciPy? >> What days are you available? So far I believe the three of us are >> available on Thursday and Friday. How about Thursday evening? Does >> anyone know what time the dinner on Thursday is? >> >> Thanks, >> >> James. >>

1 0

SciPy tutorials and talks will be recorded and posted online!
by Fernando Perez 13 Aug '09

13 Aug '09

Hi all, as you may recall, there have been recently a number of requests for videotaping the conference. I am very happy to announce that we will indeed have full video coverage this year of both tutorial tracks as well as the main talks (minus any specific talk where a speaker may object to being taped, since we'll respect such objections if they are made). Jeff Teeters and Kilian Koepsell from UC Berkeley, who have in the past made recordings like these of one of my workshops and a recent talk by Gael : http://www.archive.org/search.php?query=Fernando+Perez+scientific+python http://www.archive.org/details/ucb_py4science_2009_07_14_Gael_Varoquaux are going to perform the work. I'd like to sincerely thank: - Jeff and Kilian for offering to do this work and providing some of the recording equipment. - The Redwood Center for Theoretical Neuroscience, which provided other equipment. - Enthought, who are funding this on very short notice!!! (Especially Dave Peterson and Eric Jones, who tolerated my last minute nags very graciously). Without Enthought's last-minute support, this would simply not be happening. I really appreciate that everyone involved worked on short notice to make this possible, I hope the entire community will benefit from these resources being available. Best regards, f

1 0

Sparse complex eigensolver
by Niall Moran 12 Aug '09

12 Aug '09

Hi, I am trying to partially diagonalise a sparse complex matrix using scipy. I could not find a routine in the stable version 0.6.0 that would do it so I installed the development version 0.8.0.dev5798 and have been trying to use the routine sparse.linalg.speigs.ARPACK_eigs. It works fine when my matrix is not complex but when there are complex elements it seems to just ignore them and only take the real matrix elements. The returned eigenvalues are real (which I expect in any case) but the vectors are real also (which I know from other sources should be complex). My matrix is stored as a scipy.sparse.csr.csr_matrix (with type numpy.complex128) and I have tried passing the matvec and _mul_vector routine from this class to this function as the matrix vector multiplication routine. Any ideas or hints would be much appreciated. Many thanks, Niall.

1 0

scikits.timeseries frequencies faster than seconds?
by Alec Koumjian 12 Aug '09

12 Aug '09

I've been using the timeseries objects for a project that analyzes wind data. A lot of the data I use has a frequency of about 10 minutes, and so the package works wonderfully. However, some of the data has a frequency of 20Hz and I would like to package it inside a timeseries for the sake of uniformity. Does anyone have experience creating timeseries with these "invalid" frequencies?

2 1

Indexing question
by Michael Lerner 12 Aug '09

12 Aug '09

Hi, If I have a 1-dimensional array and I want to increment some values, I can do this: In [1]: a = array([0,0,0,0,0,0,0,0]) In [2]: b = array([1,2,4]) In [3]: a[b] += 1 In [4]: a Out[4]: array([0, 1, 1, 0, 1, 0, 0, 0]) I actually have a 3-dimensional array and a list of cells within it that I'd like to increment. I can't quite seem to get the syntax right, e.g. In [5]: a = zeros([3,3,3]) In [6]: b = array([[0,0,0],[1,0,0]]) In [7]: a[b] += 1 In [8]: a Out[8]: array([[[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]], [[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]], [[ 0., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]]]) whereas I'd like to end up with In [10]: desireda Out[10]: array([[[ 1., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]], [[ 1., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]], [[ 0., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]]]) I think speed matters too. In the worst case, I'll have an array with dimensions (400,400,30) and a list of around a million indices that I'd like to increment. Thanks, -Michael -- Michael Lerner, Ph.D. IRTA Postdoctoral Fellow Laboratory of Computational Biology NIH/NHLBI 5635 Fishers Lane, Room T909, MSC 9314 Rockville, MD 20852 (UPS/FedEx/Reality) Bethesda MD 20892-9314 (USPS)

4 4

vectorizing a function with non scalar arguments
by Luke 11 Aug '09

11 Aug '09

I have a function of the following form: def f(q, qd, parameter_list): # Unpacking the parameters m, g, I11, I22, I33 = parameter_list # Unpacking the qdots q1, q2, q3, q4, q5, q6 = q q1p, q2p, q3p, q4p, q5p, q6p = qd ....Some calculations to determine u1,u2,u3,u4,u5,u6.... return [u1, u2, u3, u4, u5, u6] As is visible above, q and qdot need to be length 6 tuples/lists, while parameter_list needs to be length 5. I would like to vectorize this function so that I could pass it: 1) q and qd as a 2-d numpy arrays of shape (n, 6) 2) parameter_list as a 1-d numpy array of shape (5,) and it would return a 2-d numpy array of shape (n, 6) parameter_list is a bunch of constants that don't change, while q and qd are the things that are different between say q[n, :] and q[n+1, :]. I realize that I can use a for loop to loop from 0 to n-1. My question is whether there a way to use vectorize (or some other decorator), to obtain the behavior? Or is there another nice method that people might be able to recommend? The documentation for vectorize is very sparse and seems like it is geared only towards function which have scalar arguments. Thanks, ~Luke

2 1

How to average different pieces or an array?
by Joseph Smidt 10 Aug '09

10 Aug '09

I'm sure this is easy I just can't think of how to do it without a bunch of for loops which I would like to avoid since they take up so much computational time. Say I have an array x, with 100 entries. Lets say I have another array y that looks like this: y = [0,5,10,20,40,100] which specified which ranges of x I would like to average. In other words, I want elements 0 - 4 in x to be replaced by their average value. 5-9 replaced with their average value. 10-19 replaced by their averaged value, etc... Effectively this will give me a new array with one 100 entries with different sized binning along the array. If anyone knows how to do this without using for loops the whole way I would appreciate it. Example: I will make it simpler in case my above description was confusing. Say x = [1,3,2,6,7,4,5,4,9,4] and I passed in y = [0,2,4,10]. The above description would return a new X: X = [2,2,4,4,5.5,5.5,5.5,5.5,5.5,5.5] So that the apprpriate bins are averaged. Joseph Smidt -- ------------------------------------------------------------------------ Joseph Smidt <josephsmidt(a)gmail.com> Physics and Astronomy 4129 Frederick Reines Hall Irvine, CA 92697-4575 Office: 949-824-3269

6 10