[SciPy-Dev] Find points in delaunay triangulation : scipy.spatial vs. scipy.interpolation

[Pablo Winant]

Hi,

I tried to use interpolation routines in scipy recently and I have found 
two slight performance issues

  - The LinearNDInterpolation object implemented in cython requires a 
list of points and a list of values to be created. But is is not 
documented how to change the values of the interpolator without doing 
the mesh again. This is useful when one is solving the values of a 
function at the vertices of the mesh : one doesn't want to do the 
triangulation again and again. Maybe there could be a simple specific 
method to set the values in this case. In that case it would consist in 
changing the value of a property but it would be consistent with more 
general interpolation schemes.

- I tried to use the delaunay object from scipy and noticed a strange 
thing: for a given set of coordinates it takes longer to get the indices 
of the triangles containing the points than it takes to perform the 
interpolation using LinearND object. This is puzzling since apparently 
the implementation of LinearND performs many calls to the qhull library 
to get this indices. Attached is a simple exampe demonstrating this anomaly.

One last thing: I have written an interpolation object on sparse grids, 
using smolyak product of chebychev polynomials. It is written in pure 
python (vectorized) and licensed under the bsd license. Currently it 
lives in another library but I guess it would make more sense to have 
something like that in a more general scientific lib. Let me know if you 
are interested. (it is available there anyway: 
https://github.com/albop/dynare-python/tree/master/dolo/src/dolo/numeric: chebychev.py 
and smolyak.py)

Best regards,

Pablo