Mailman 3 August 2009 - SciPy-User

Installing latest svn version UMFPACK issue
by Pim Schellart 04 Aug '09

04 Aug '09

Dear Scipy users/developers, I want to install the latest svn version of scipy because I need the curve_fit routine in minpack. This routine (although not very robust) makes curve fitting much easier, and I would like to see it included in the release version of scipy soon :) However when I try to install the latest svn checkout I get an error that the installation cannot find UMFPACK on my system (OSX). I did not get this error with previous installations and I am wondering why I do get it now. Is this a newly added dependency (in which case please tell me how to install it since it seems to be designed to be integrated with matlab), if not why does this problem occur (is there perhaps something wrong with the installation script). Kind regards, Pim Schellart

3 2

Dimensions of x and y are incompatible
by Omer Khalid 04 Aug '09

04 Aug '09

Hi, I am trying to plot two curves on a same graph. The two input datasets are xdata, ydata and they range between (0.01, 1] for xdata and (0.1, 1] for ydata. I am constantly getting this error when plotting "Dimensions of x and y are incompatible" plt.figure() xdata= array(xd) ydata = array(xd) N=len(xdata) print "N: %d" % (N) x=arange(0.0, 1.0+0.01, 0.01) M=len(ydata) print "M: %d" % (M) y=arange(0.0, 1.0+0.1, 0.1) plt.plot(x,xdata,'r') plt.plot(y,ydata,'b') plot.legend((xl,yl)) plt.grid(True) plt.title(title) plt.xlabel('Value Ratio') plt.ylabel('Simulation Time Unit') pdf_name = '%s-%s.pdf' % (str(filename),strftime("%m%d-%M-%S")) pwd = os.getcwd() path = os.path.join(pwd, exp, pdf_name) plt.savefig(path) I have also tried: x=arange(1,N) y=arange(1,M) But it didn't work. Any help would be much appreciated. Thanks, Omer

4 3

optical ray tracing
by Elliot Hallmark 04 Aug '09

04 Aug '09

howdy, I'm looking to do some raytracing for a nonimaging design and there are no opensource packages that address the geometry needed. Scipy seems like a good package to become familiar with anyway (I miss mathematica now that i'm no longer a university student), but i was wondering if anyone has already done some optical ray tracing work in scipy that i could build off of. Wouldn't be to hard to do from scratch, but my focus is on the design and the fewer excursions into building new tools the better. thanks, elliot

7 8

Scipy and statistics: probability density function
by Daniel J Farrell 04 Aug '09

04 Aug '09

Dear list, I am looking for some of the functionality provided by the GNU Scientific Library histogram module (http://www.gnu.org/software/gsl/manual/html_node/Histograms.html ). In particular, a need to be able to create a probability density function from my histogram of data. This will allow inverse look-ups to be performed, i.e. for a random number (0-->1) find the associated probability (http://www.gnu.org/software/gsl/manual/html_node/The-histogram-probability-… ). This allows the distribution and for samples to be returned weighted by the probability of the distribution -- which is a common task! Is is easier just to call GSL using C-types? I just wanted to check that this functionality doesn't already exist before I stray too far, or implement it my self. Best regards Dan

3 3

expectation function for stats.distributions
by josef.pktd＠gmail.com 03 Aug '09

03 Aug '09

I just needed to verify an expectation for fixing of the stats.models code. A while ago a wrote a function that can be attached as method to the distributions classes. It's pretty simple, but comes in handy when an expectation or conditional expectation has to be checked. Is there an interest in adding this as a new method to the distributions? explanations are here http://jpktd.blogspot.com/2009/04/having-fun-with-expectations.html copy of file below Josef ''' copy from webpage, where is the original? Author jpktd ''' import numpy as np from scipy import stats, integrate def expectedfunc(self, fn=None, args=(), lb=None, ub=None, conditional=False): '''calculate expected value of a function with respect to the distribution only for standard version of distribution, location and scale not tested Parameters ---------- all parameters are keyword parameters fn : function (default: identity mapping) Function for which integral is calculated. Takes only one argument. args : tuple argument (parameters) of the distribution lb, ub : numbers lower and upper bound for integration, default is set to the support of the distribution conditional : boolean (False) If true then the integral is corrected by the conditional probability of the integration interval. The return value is the expectation of the function, conditional on being in the given interval. Returns ------- expected value : float ''' if fn is None: def fun(x, *args): return x*self.pdf(x, *args) else: def fun(x, *args): return fn(x)*self.pdf(x, *args) if lb is None: lb = self.a if ub is None: ub = self.b if conditional: invfac = self.sf(lb,*args) - self.sf(ub,*args) else: invfac = 1.0 return integrate.quad(fun, lb, ub, args=args)[0]/invfac print expectedfunc(stats.norm, lambda(x): (x)**4) print expectedfunc(stats.norm, lambda(x): min((x)**2,1.5**2)) #Jonathans version from scipy.stats import norm as Gaussian c = 1.5 # our "cutoff" point c = 0.5 # try another value tmp = 2 * Gaussian.cdf(c) - 1 gamma = tmp + c**2 * (1 - tmp) - 2 * c * Gaussian.pdf(c) print gamma print expectedfunc(stats.norm, lambda(x): min((x)**2,c**2))

3 5

An extra parameter to stats.chisquare ?
by Pierre GM 03 Aug '09

03 Aug '09

All, stats.chisquare requires a mandatory parameter (frequency of observations) and an optional argument (theoretical frequencies). In that second case, I think we have to introduce yet another parameter, p, corresponding to the number of parameters of the theoretical distribution: the number of degrees-of-freedom for the chisqprob would then be k-p (with k the sample size), instead of the current k-1. Of course, we can set p=1 by default. Comments ? P.

2 6

better cumulative integration routines (like cumtrapz)
by Roban Hultman Kramer 03 Aug '09

03 Aug '09

Hi, I'm wondering if anyone can recommend a routine for one-dimensional cumulative numerical integrals. Obviously there is scipy.integrate.cumtrapz, but I was wondering if anyone knows of an implementation of a more sophisticated numerical integration algorithm that can take not just end points, but a whole series of points, and efficiently approximate the definite integral at each point. In other words, I have a function f(x), and want the value of the definite integral from x0 to [x1, x2, x3, x4, ...] Ideally I would want to be able to specify the desired precision of the results, and also get an estimate of the achieved precision. Essentially I'm asking for scipy.integrate.quad, but optimized for a whole set of points. For the most part my functions are well behaved, so I don't need fancy singularity-handling, though built-in infinite-limit handling is great. The obvious hack is to compute the definite integral separately for each interval in my series of points, then add them together, but I was wondering if anyone knows of an algorithm/implementation that does this efficiently. I'd be willing to work on implementing something myself, if you have pointers to an algorithm. Thanks! -Roban

1 0

SciPy Foundation
by Joe Harrington 03 Aug '09

03 Aug '09

About sixteen months ago, I launched the SciPy Documentation Project and its Marathon. Dozens pitched in and now numpy docs are rapidly approaching a professional level. The "pink wave" ("Needs Review" status) is at 56% today! There is consensus among doc writers that much of the rest can be labeled in the "unimportant" category, so we're close to starting the review push (hold your fire, there is a web site mod to be done first). We're also nearing the end of the summer, and it's time to look ahead. The path for docs is clear, but the path for SciPy is not. I think our weakest area right now is organization of the project. There is no consensus-based plan for improvement of the whole toward a stated goal, no centralized coordination of work, and no funded work focused on many of our weaknesses, notwithstanding my doc effort and what Enthought does for code. I define success as popular adoption in preference to commercial packages. I believe in vote-with-your-feet: this goal will not be reached until all aspects of the package and its presentation to the world exceed those of our commercial competition. Scipy is now a grass roots effort, but that takes it only so far. Other projects, such as OpenOffice and Sage, don't follow this model and do produce quality products that compete with commercial offerings, at least on open-source platforms. Before we can even hope for that, we have to do the following: - Docs - Rest of numpy reference pages reviewed and proofed or marked unimportant - Scipy reference pages - User manual for the whole toolstack - Multiple commercial books - Packaging - Personal Package Archive or equivalent for every release of every OS for the full toolstack (There are tools that do this but we don't use them. NSF requires Metronome - http://nmi.cs.wisc.edu/ - for funding most development grants, so right now we're not even on NSF's radar.) - Track record of having the whole toolstack installation "just work" in a few command lines or clicks for *everyone* - Regular, scheduled releases of numpy and scipy - Coordinated releases of numpy, scipy, and stable scikits into PPA system - Public communication - A real marketing plan - Executing on that plan - Web site geared toward multiple audiences, run by experts at that kind of communication - More webinars, conference booths, training, aimed at all levels - Demos, testimonials, topical forums, all showcased - Code - A full design review for numpy 2.0 - No more inconsistencies like median(), lacking "out", degrees option for angle functions? - Trimming of financial functions, maybe others, from numpy? - Package structure review (eliminate "fromnumeric"?) - Goal that this be the last breakage for numpy API (the real 1.0) - Scipy - Is it maintainable? should it be broken up? - Clear code addition path (or decide never to add more) - Docs (see above) - Add-on packages - Both existence of and good indexing/integration/support for field-specific packages - Clearer development path for new packages - Central hosting system for packages (svn, mailing lists, web, build integration, etc.) - Simultaneous releases of stable packages along with numpy/scipy I posted a basic improvement plan some years back. The core ideas have not changed; it is linked from the bottom of http://scipy.org/Developer_Zone. I chose our major weakness to begin with and started the doc project, using some money I could justify spending simply for the utility of docs for my own research. I funded the work of two doc coordinators, one each this summer and last. Looking at http://docs.scipy.org/numpy/stats/, you can see that when a doc coordinator was being paid (summers), work got done. When not, then not. Without publicly announcing what these guys made, I'll be the first to admit that it wasn't a lot. Yet, those small sums bought a huge contribution to numpy through the work of several dozen volunteers and the major contributions of a few. My conclusion is that active and constant coordination is central to motivating volunteer work, and that without a salary we cannot depend on coordination remaining active. On the other hand, I have heard Enthought's leaders bemoan the high cost of devoting employee time to this project, and the low returns available from selling support to universities and non-profit research institutes. Their leadership has moved us forward, particularly in the area of code, but has not provided the momentum necessary to carry us forward on all fronts. It is time for the public and education sectors to kick in some resources and organizational leadership. We are, after all, benefitting immensely. Since the cost of employee time is not so high for us in the public and education sectors, I propose to continue hiring people like Stefan and David as UCF employees or contractors, and to expand to hiring others in areas like packaging and marketing, provided that funding for those hires can be found. However, my grant situation is no longer as rich as it has been the past two years, and the needs going forward are greater than in the past if we're now to tackle all the points above. So, I will not be hiring another doc guru from my research grants next year. I am confident that others are willing to pitch in financially, but few will pitch in a full FTE, and we need several. We can (and will) set up a donations site, but donation sites tend to receive pizza money unless a sugar daddy comes along. Those benefitting most from the software, notably education, non-profit research, and government institutions, are *forbidden* from making donations by the terms of their grants. NSF doesn't give you money so you can give it away. We need to provide services they can buy on subcontract and a means for handling payments from them. Selling support does not solve the problem, as that requires spending most of the income on servicing that particular client. Rather, we need to sell a chunk of documentation or the packaging of a particular release, and then provide the product not just to that client but to everyone. We can also propose directly for federal and corporate grant funds. I have spoken with several NASA and NSF program managers and with Google's Federal Accounts Representative, and the possibilities for funding are good. But, I am not going to do this alone. We need a strong proposal team to be credible. So, I am seeking a group that is willing to work with me to put up the infrastructure of a funded project, to write grant proposals, and to coordinate a financial effort. Members of this group must have a track record of funded grants, business success, foundation support, etc. We might call it the SciPy Foundation. It could be based at UCF, which has a low overhead rate and has infrastructure (like an HR staff), or it might be independent if we can find a good director willing to devote significant time for relatively low pay compared to what they can likely make elsewhere. I would envision hiring permanent coordinators for docs, packaging, and marketing communications. Enthought appears to have code covered by virtue of having hired Travis, Robert, etc.; how to integrate that with this effort is an open question but not a difficult one, I think, as code is our strongest asset at this point. I invite discussion of this approach and the task list above on the scipy-dev(a)scipy.org mailing list. If you are seeing this post elsewhere, please reply only on scipy-dev(a)scipy.org. If you are eligible to lead funding proposals and are interested in participating in grant writing and management activities related to work in our weak areas, please contact me directly. Thanks, --jh-- Prof. Joseph Harrington Planetary Sciences Group Department of Physics MAP 414 4000 Central Florida Blvd. University of Central Florida Orlando, FL 32816-2385 jh(a)physics.ucf.edu planets.ucf.edu

5 9

works in gnuplot, but not quite in Python/Gnuplot.py - splot (sin(x)/x) * (sin(y)/y)
by Esmail 03 Aug '09

03 Aug '09

Hello all, I am trying to generate a plot with Python/Gnuplot.py that I can create interactively with gnuplot. ---------- gnuplot> set isosample 40 gnuplot> set hidden3d gnuplot> splot (sin(x)/x) * (sin(y)/y) ---------- gives me a *beautiful* continuous plot. I am trying to do the same with Python/Gnuplot.py - but no joy. I have the following short Python code to generate my data which I then funnel into a file called 'somb.dat' ---------- #!/usr/bin/env python # somb.py # compute the sombrero function from math import sin from numpy import arange def f(x, y): return (sin(x)/x) * (sin(y)/y) for x in arange(-10, 11): for y in arange(-10, 11): print '%3d %3d %f' % (x, y, f(x, y)) print ---------- e.g., somb.py > somb.dat and I use this to plot the data: ---------- #!/usr/bin/env python import Gnuplot g = Gnuplot.Gnuplot(debug=1) g('set hidden3d') g('set isosample 40') g('set datafile missing "nan"') g.splot(Gnuplot.File('somb.dat', with_='lines linewidth 1')) raw_input() # to pause ---------- Is there a way to achieve the same sort of continuous output that gnuplot generates? I tried using 0.1 as 3rd parameter for arange, but it didn't help. I suspect the problem is with the nan values (and/or the numeric representation in Python)?. How does gnuplot handle this? How can I achieve the same with Python/Gnuplot.py? Any suggestions/hints? (I am thrilled that I can use gnuplot via Python) Thanks, Esmail

1 1

sorting an array
by Martin 03 Aug '09

03 Aug '09

Hi, I would ideally like to sort an multidimensional array by the first column, but I can only manage to adjust the whole array. For example... >a = array([(1,0,2,5),(1,5,6,2),(1,2,8,7),(1,4,4,6),(1,3,2,3),(1,11,2,0),(1,10,1,3),(1,9,0,4),(1,8,9,6)]) >a array([[ 1, 0, 2, 5], [ 1, 5, 6, 2], [ 1, 2, 8, 7], [ 1, 4, 4, 6], [ 1, 3, 2, 3], [ 1, 11, 2, 0], [ 1, 10, 1, 3], [ 1, 9, 0, 4], [ 1, 8, 9, 6]]) >sort(a, axis = 0) array([[ 1, 0, 0, 0], [ 1, 2, 1, 2], [ 1, 3, 2, 3], [ 1, 4, 2, 3], [ 1, 5, 2, 4], [ 1, 8, 4, 5], [ 1, 9, 6, 6], [ 1, 10, 8, 6], [ 1, 11, 9, 7]]) When really what i wanted was not for the 3rd and 4th columns to also be sorted nunmerically in that way. In unix I would reformat the first two columns in gawk i.e. gawk '{printf("%2.3d %2.3d, %d, %d\n", $1, $2, $3, $4)}' | sort -n -k 1 Any help would be appreciated! Thanks Martin

5 10