Mailman 3 September 2010 - SciPy-Dev

Unexpected behaviour from timeseries.filled (in Python 2.7)
by Åsmund Hjulstad Sept. 6, 2010

Sept. 6, 2010

It seems something in Python 2.7 / Numpy 1.5 (or my Python 2.7 installation) is breaking scikits.timeseries.filled import scikits.timeseries from numpy import array t3 = scikits.timeseries.time_series(array([0.]), mask=[True],dates=['2010-Aug'],freq='M') t3.fill_value= 0. print(t3) print(t3.data) print(t3.fill_value) print(t3.filled()) # It seems it ignores the fill_value parameter print(t3.filled(0.)) In Python 2.6 this outputs c:\python26\python scikits_bug.py [--] [ 0.] 0.0 [ 0.] [ 0.] , but in Python 2.7 I get: c:\python27\python scikits_bug.py [--] [ 0.] 0.0 [ 1.00000000e+20] # This must be wrong! [ 0.] Python 2.6: scikits.timeseries version 0.91.3, numpy 1.3.0 Python 2.7: scikits.timeseries version 0.91.3, numpy 1.5.0 -- Åsmund Hjulstad, asmund(a)hjulstad.com

2 4

N-dimensional interpolation
by Pauli Virtanen Sept. 6, 2010

Sept. 6, 2010

Hi all, I took the Qhull by the horns, and wrote a straightforward `griddata` implementation for working in N-D: http://github.com/pv/scipy-work/tree/qhull http://github.com/pv/scipy-work/blob/qhull/scipy/interpolate/qhull.pyx Sample (4-D): ------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import griddata def f(x, y, z, u): return np.cos(x + np.cos(y + np.cos(z + np.cos(u)))) points = np.random.rand(2000, 4) values = f(*points.T) p = np.linspace(0, 1, 1500) z = griddata(points, values, np.c_[p, p, p, p], method='linear') plt.plot(p, z, '-', p, f(p, p, p, p), '-') plt.show() ------------------------------------------------------------- It performs the N-D Delaunay tesselation with Qhull, and after that does linear barycentric interpolation on the simplex containing each point. (The nearest-neighbor "interpolation" is, on the other hand, implemented using scipy.spatial.KDTree.) I ended up writing some custom simplex-walking code for locating the simplex containing the points -- this is much faster than a brute-force search. However, it is slightly different from what `qh_findbestfacet` does. [If you're a computational geometry expert, see below...] Speed-wise, Qhull appears to lose in 2-D to `scikits.delaunay` by a factor of about 10 in triangulation in my tests; this, however, includes the time taken by computing the barycentric coordinate transforms. Interpolation, however, seems to be faster. *** I think this would be a nice addition to `scipy.optimize`. I'd like to get it in for 0.9. Later on, we can specialize the `griddata` function to perform better in low dimensions (ndim=2, 3) and add more interpolation methods there; for example natural neighbors, interpolating splines, etc. In ndim > 3, the `griddata` is however already feature-equivalent to Octave and MATLAB. Comments? Pauli PS. A question for computational geometry experts: `qh_findbestfacet` finds the facet on the lower convex hull whose oriented hyperplane distance to the point lifted on the paraboloid is maximal --- however, this does not always give the simplex containing the point in the projected Delaunay tesselation. There's a counterexample in qhull.pyx:_find_simplex docstring. On the other hand, Qhull documentation claims that this is the way to locate the Delaunay simplex containing a point. What gives? It's clear, however, that if a simplex contains the point, then the hyperplane distance is positive. So currently, I just walk around positive-distance simplices checking the inside condition for each and hoping for the best. If nothing is found, I just fall back to brute force search. This seems to work well in practice. However, if a point is outside the triangulation (but inside the bounding box), I have to fall back to a brute-force search. Is there a better option here? *** Another sample (2-D, object-oriented interface): ------------------------------------------------------------- import time import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import LinearNDInterpolator # some random data x = np.array([(0,0), (0,1), (1, 1), (1, 0)], dtype=np.double) np.random.seed(0) xr = np.random.rand(200*200, 2).astype(np.double) x = np.r_[xr, x] y = np.arange(x.shape[0], dtype=np.double) # evaluate on a grid nn = 500 xx = np.linspace(-0.1, 1.1, nn) yy = np.linspace(-0.1, 1.1, nn) xx, yy = np.broadcast_arrays(xx[:,None], yy[None,:]) xi = np.array([xx, yy]).transpose(1,2,0).copy() start = time.time() ip = LinearNDInterpolator(x, y) print "Triangulation", time.time() - start start = time.time() zi = ip(xi) print "Interpolation", time.time() - start from scikits.delaunay import Triangulation start = time.time() tri2 = Triangulation(x[:,0], x[:,1]) print "scikits.delaunay triangulation", time.time() - start start = time.time() ip2 = tri2.linear_interpolator(y) zi2 = ip2[-0.1:1.1:500j,-0.1:1.1:500j].T print "scikits.delaunay interpolation", time.time() - start print "rel-difference", np.nanmax(abs(zi - zi2))/np.nanmax(np.abs(zi)) plt.figure() plt.imshow(zi) plt.clim(y.min(), y.max()) plt.colorbar() plt.figure() plt.imshow(zi2) plt.clim(y.min(), y.max()) plt.colorbar() plt.show() ------------------------------------------------------------- -- Pauli Virtanen

5 8

Nonlinear large-scale root finding
by Pauli Virtanen Sept. 5, 2010

Sept. 5, 2010

Hi, I finally committed the scipy.optimize.nonlin rewrite that was pending for a long time. What is now there for large-scale root finding is: * Newton-Krylov method (``scipy.optimize.newton_krylov``) * Generalized secant methods: - Limited-memory Broyden (``scipy.optimize.broyden1``, ``scipy.optimize.broyden2``) - Anderson (``scipy.optimize.anderson``) * Simple iterations (``scipy.optimize.diagbroyden``, ``scipy.optimize.excitingmixing``, ``scipy.optimize.linearmixing``) Here's a simple demo: http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#large-prob… -- Pauli Virtanen

1 0

Re: [SciPy-Dev] RFR: N-dimensional interpolation
by Eric Bruning Sept. 5, 2010

Sept. 5, 2010

I made a trial run with some 3D data from a weather radar this evening. I'm interpolating to a horizontal slice through a (distorted) spherical grid, as shown here: http://bitbucket.org/deeplycloudy/pyrsl/src/tip/example.py Both linear and nearest-neighbor work great. Image of linear interpolation at http://yfrog.com/msspop Many thanks for cooking this up. I've been hoping for something like this for a long time. -Eric

1 0

Scipy.org down (again)
by Fernando Perez Sept. 1, 2010

Sept. 1, 2010

Howdy, the usual... Cheers, f

1 0