Mailman 3 March 2008 - SciPy-User

NumPy (1.0.5) DocDay (Fri, Mar. 21)
by Jarrod Millman 20 Mar '08

20 Mar '08

Hello, As I mentioned yesterday, I am holding a NumPy DocDay on Friday, March 21st. I am in Paris near the RER B or C Saint-Michel station (with Stefan van der Walt, Matthieu Brucher, and Gael Varoquaux). If you are in the area and want to join us just send me an email by the end of tonight and I will let you know where we are meeting. If you can't stop by, but are still willing to help out we will convene on IRC during the day on Friday (9:30am-?? GMT+1). Come join us at irc.freenode.net (channel scipy). We may update the list of priorities which is still located on the NumPy Trac Wiki: http://projects.scipy.org/scipy/numpy/wiki/DocDays While I am hoping to have everyone focus on NumPy, I would be happy if anyone wants to work on SciPy documentation as well: http://projects.scipy.org/scipy/scipy/wiki/DocDays Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

1 0

Splines for even functions?
by Michael McNeil Forbes 20 Mar '08

20 Mar '08

Does anyone have a suggestion on how to force a smoothing spline to be an even function f(x) = f(-x)? I am trying to interpolate the following data representing an even function and I need a smooth but even spline. The problem is that, in the process of smoothing, the knots are not chosen evenly so the resulting spline is not even. Increasing the tolerances ultimately forces more knots to be added resulting in an even spline, but it is difficult to automate this process. from scipy import linspace, interpolate, array X = array([-1. , -0.65016502, -0.58856235, -0.26903553, -0.17370892, -0.10011001, 0. , 0.10011001, 0.17370892, 0.26903553, 0.58856235, 0.65016502, 1.]) Y = array([ 1. , 0.62928599, 0.5797223 , 0.39965815, 0.36322694, 0.3508061 , 0.35214793, 0.3508061 , 0.36322694, 0.39965815, 0.5797223 , 0.62928599, 1.]) E = array([ 1.00000000e-12, 1.45164012e-03, 2.04367440e-03, 2.34266209e-03, 1.64542215e-03, 2.21561749e-03, 3.14980262e-03, 2.21561749e-03, 1.64542215e-03, 2.34266209e-03, 2.04367440e-03, 1.45164012e-03, 1.00000000e-12]) f = interpolate.UnivariateSpline(X,Y,1./E,k=3) x = linspace(-1,1,10) max(abs(f(x) - f(-x))) 0.00066571302580914482 Thanks, Michael.

1 1

JOB: Short-term programming (consultant) work
by Rob Clewley 19 Mar '08

19 Mar '08

Dear SciPy users, The developers of the PyDSTool dynamical systems software project have money to hire a Python programmer on a short-term, per-task basis as a technical consultant. The work can be done remotely and will be paid after the completion of project milestones. The work must be completed by July, when the current funds expire. Prospective consultants could be professionals or students and will have proven experience and interest in working with NumPy/SciPy, scientific computation in general, and interfacing Python with C and Fortran codes. Detailed work plan, schedule, and project specs are negotiable (if you are talented and experienced we would like your input). The rate of pay is commensurate with experience, and may be up to $45/hr or $1000 per project milestone (no fringe benefits), according to an agreed measure of satisfactory product performance. There is a strong possibility of longer term work depending on progress and funding availability. PyDSTool (pydstool.sourceforge.net) is a multi-platform, open-source environment offering a range of library tools and utils for research in dynamical systems modeling for scientists and engineers. As a research project, it presently contains prototype code that we would like to improve and better integrate into our long-term vision and with other emerging (open-source) software tools. Depending on interest and experience, current projects might include: * Conversion and "pythonification" of old Matlab code for model analysis * Improved interface for legacy C and Fortran code (numerical integrators) via some combination of SWIG, Scons, automake * Overhaul of support for symbolic processing (probably by an interface to SymPy) For more details please contact Dr. Rob Clewley (rclewley) at (@) the Department of Mathematics, Georgia State University (gsu.edu). -- Robert H. Clewley, Ph. D. Assistant Professor Department of Mathematics and Statistics Georgia State University 720 COE, 30 Pryor St Atlanta, GA 30303, USA tel: 404-413-6420 fax: 404-651-2246 http://www.mathstat.gsu.edu/~matrhc http://brainsbehavior.gsu.edu/

1 0

Reference for scipy.stats.kurtosistest
by Kurt Smith 19 Mar '08

19 Mar '08

Hello scipy-users (Specifically Robert Kern): I'm making extensive use of scipy.stats.kurtosis, and noticed 'kurtosistest' which calculates the Z-score, and 2-tail Z-probability. The source for this function is opaque to me -- what is being calculated, and how? Is there a reference to which you can refer me? Thanks for your time, Kurt

2 2

Incremental histogram?
by Roger Herikstad 19 Mar '08

19 Mar '08

Hi list, I need to histogram an array of long ints, but the array itself is too big to keep in memory. I was thinking of using an incremental approach, i.e. assign each sample in the array to the appropriate bin, sample by sample. Right now, I have the array (well, list really) constructed as a generator, and I was wondering if anyone has an efficient algorithm for doing histogram count on such a generator object? ~ Roger

6 8

Re: [SciPy-user] Runge-Kutta ODE vs. built-in tools
by Warren Weckesser 19 Mar '08

19 Mar '08

Zane (and everyone else), You might be interested in trying VFGEN (www.warrenweckesser.net/vfgen). I just released version 2.2.0. VFGEN is a program that takes an XML description of a vector field (i.e. differential equations) and generates code for a wide variety of numerical software packages, including SciPy. A simple example is the damped pendulum: θ' = v v' = -b*v/(m*L^2) - (g/L)*sin(θ) An XML vector field file for this system is ---------------------------------------------------------------- <?xml version="1.0" ?> <VectorField Name="pendulum" Description="Pendulum Vector Field"> <Parameter Name="g" Description="gravitational constant" DefaultValue="9.81" /> <Parameter Name="b" Description="friction constant" DefaultValue="0.0" /> <Parameter Name="L" Description="pendulum length" DefaultValue="1.0" /> <Parameter Name="m" Description="mass" DefaultValue="1.0" /> <StateVariable Name="theta" Description="Angle, measured from straight down" Formula="v" PeriodFrom="0" PeriodTo="2*Pi" DefaultInitialCondition="Pi-0.01" /> <StateVariable Name="v" Description="angular velocity" Formula="-b*v/(m*L^2)-(g/L)*sin(theta)" DefaultInitialCondition="0.0" /> <Function Name="energy" Description="total energy (kinetic plus potential)" Formula="m*L^2*v^2/2 - m*g*L*cos(theta)" /> </VectorField> ---------------------------------------------------------------- The VFGEN scipy command reads this file and generates this Python code: ---------------------------------------------------------------- # # pendulum.py # # Python file for the vector field named: pendulum # This file implements the vector field as # the function pendulum_vf. This function # can be used with the SciPy ODEINT function. # # This file was generated by the program VFGEN (Version:2.2.0) # Generated on 17-Mar-2008 at 00:48 # from math import * import numpy # # The vector field. # def pendulum_vf(y_,t_,p_): """ The vector field function for the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] f_ = numpy.zeros((2,)) f_[0] = v f_[1] = -sin(theta)/L*g-1.0/m/(L*L)*b*v return f_ # # The Jacobian. # def pendulum_jac(y_, t_, p_): """ The Jacobian of the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] # Create the Jacobian matrix: jac_ = numpy.zeros((2,2)) jac_[0,1] = 1.0 jac_[1,0] = -1.0/L*cos(theta)*g jac_[1,1] = -1.0/m/(L*L)*b return jac_ # # User function: energy # def pendulum_energy(y_, t_, p_): """ The user-defined function "energy" for the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] return m*(L*L)*(v*v)/2.0-m*L*cos(theta)*g ---------------------------------------------------------------- VFGEN can also generate a simple command-line solver written in Python that uses the above function to solve the initial value problem. Take a look here: http://www.warrenweckesser.net/vfgen/menu_scipy.html If you decide you would like to use a different package (e.g. MATLAB, Scilab, C with CVODE or with the GNU Scientific Library, Fortran with LSODA or RADAU5, etc), another one-line VFGEN command will generate the code that you need. Best regards, Warren Weckesser

1 0

extending a 2-dim array ?
by Stef Mientki 18 Mar '08

18 Mar '08

hello, To combine a number of signals (of possible different length), I need to extend a 2-dimensional array. I succeeded, but it looks much too complex in my opinion: a_list = [1,2,3,4,5] a_vector = array ( a_list ) b_list = [11,22,33,44,55] b_vector = array ( b_list ) ab_array = vstack (( a_vector, b_vector )) # extend each signal with 1 zero element extend = transpose ( vstack ( ( transpose ( ab_array ), zeros ((1,2)) )) ) does anyone know a more elegant (= more simple) method ? thanks, Stef Mientki

3 2

normalized frequencies to pi radians / sample
by Forrest Sheng Bao 18 Mar '08

18 Mar '08

Hi, I am trying to use buttord function to compute the order and frequency of a Butterworth filter. But I got confused about the passband and stopband edge frequency. According to the doc, wp, ws -- Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample). For example: Lowpass: wp = 0.2, ws = 0.3 Highpass: wp = 0.3, ws = 0.2 Bandpass: wp = [0.2, 0.5], ws = [0.1, 0.6] Bandstop: wp = [0.1, 0.6], ws = [0.2, 0.5] I should normalize the frequency to pi radians / sample. Now, suppose a frequency 10Hz and my sampling rate is 100 Hz. Is the normalized frequency 0.1? Thanks, Forrest -- Forrest Sheng Bao Ph.D. student, Dept. of Computer Science M.Sc. student, Dept. of Electrical & Computer Engineering Rm 115, Experimental Sciences Building Texas Tech University, USA http://fsbao.net 1-806-577-4592 Forrest is an equal opportunity Email sender. 1. You are encouraged to use the language you prefer. Beyond English, I can also read traditional/simplified Chinese and a bit German. 2. I will only send you files readable to free or open source software.

3 2

timeseries and maskedarrays
by bevan 18 Mar '08

18 Mar '08

Hello, I recently posted on the numpy forum looking for a solution to creating totals and averages for time series of rainfall. It was suggested to me that I try the timeseries package. I have managed to get it working (struggling due to limited abilities rather than any issues with the package) - thanks by the way to the developers. However, when I run my code i get the following warning: C:\Python25\lib\site-packages\scipy\sandbox\maskedarray\core.py:1521: UserWarning: Warning: converting a masked element to nan. warnings.warn("Warning: converting a masked element to nan.") Is this an issue? The main question I have is: How would i create an average on values based on month, that is sum the rainfall from daily values to monthly totals (done), then average all the Januarys, Februarys etc in the timeseries? Also any other issues with my code (below) that is likely to trip me up later? Thanks for your time. import time import numpy import maskedarray as MA import timeseries as TS def readTSFtimeseries(filename): """ # PURPOSE: read a HYDROL .tsf file into a timeseries array. # # INPUT: #20/11/2003 @ 00:00:00 Q255 T5 #21/11/2003 @ 00:00:00 0.000000 Q30 T5 #............................................ #07/05/2008 @ 00:00:00 2.500000 Q10 T5 # # OUTPUT: a timeseries array # # EXAMPLE: >>> data = readTSFtimeseries(r'B:\Python\test.tsf') # >>> x = data.dates # >>> y = data.series """ MYFILE = open(filename, "r") datafile = MYFILE.readlines() MYFILE.close() DateList=[] Rain_mmList=[] QualityList=[] TypeList=[] jd=[] for datalines in datafile: dataflds = datalines.split() datefld =time.strptime(dataflds[0]+'_'+dataflds[2], '%d/%m/%Y_%H:%M:%S') jd.append(datefld[7]) DateList.append(TS.Date('D',year=datefld[0],month=datefld[1],day=datefld [2])) if len(dataflds) == 6: Rain_mmList.append(float(dataflds[3])) QualityList.append(dataflds[4]) TypeList.append(dataflds[5]) if len(dataflds) == 5: #Rain_mmList.append(float(-9.99)) Rain_mmList.append(TS.tsmasked) QualityList.append(dataflds[3]) TypeList.append(dataflds[4]) DateArr=TS.DateArray(dates=DateList,freq='D',copy =False) data=TS.time_series(Rain_mmList,DateArr) print data.dates.size,data.series.size return(data) Raindata =readTSFtimeseries(r"C:\Documents and Settings\bevanj\Desktop\rain.tsf") import pylab import matplotlib from timeseries import plotlib as TSPL MonRaindata = Raindata.convert('monthly', MA.sum) print MonRaindata print Raindata.convert('yearly', MA.sum) #fig1=TSPL.tsfigure()3 #fplt1=fig1.add_tsplot(111) #fplt1.tsplot(RainData,'-') #fplt1.tsplot(MonRainData,'-') #pylab.show()

4 4

ANN: EuroSciPy 2008 Conference - Leipzig, Germany
by Travis Vaught 17 Mar '08

17 Mar '08

Greetings, We're pleased to announce the EuroSciPy 2008 Conference to be held in Leipzig, Germany on July 26-27, 2008. http://www.scipy.org/EuroSciPy2008 We are very excited to create a venue for the European community of users of the Python programming language in science. This conference will bring the presentations and collaboration that we've enjoyed at Caltech each year closer to home for many users of SciPy, NumPy and Python generally--with a similar focus and schedule. Call for Participation: ---------------------- If you are a scientist using Python for your computational work, we'd love to have you formally present your results, methods or experiences. To apply to present a talk at this year's EuroSciPy, please submit an abstract of your talk as a PDF, MS Word or plain text file to euroabstracts(a)scipy.org. The deadline for abstract submission is April 30, 2008. Papers and/or presentation slides are acceptable and are due by June 15, 2008. Presentations will be allotted 30 minutes. Registration: ------------ Registration will open April 1, 2008. The registration fee will be 100.00€ for early registrants and will increase to 150.00€ for late registration. Registration will include breakfast, snacks and lunch for Saturday and Sunday. Volunteers Welcome: ------------------ If you're interested in volunteering to help organize things, please email us at info(a)scipy.org.

9 18