Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

alternate fix for issue #997 #1477

Merged
merged 6 commits into from Dec 4, 2012
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
4 changes: 2 additions & 2 deletions lib/matplotlib/delaunay/_delaunay.cpp
Expand Up @@ -289,8 +289,8 @@ static PyObject *linear_interpolate_grid(double x0, double x1, int xsteps,
if (!z) return NULL;
z_ptr = (double*)PyArray_DATA(z);

dx = (x1 - x0) / (xsteps-1);
dy = (y1 - y0) / (ysteps-1);
dx = ( x1==x0 ? 0 : (x1 - x0) / (xsteps-1) );
dy = ( y1==y0 ? 0 : (y1 - y0) / (ysteps-1) );
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps I am missing something here, but if x1==x0, then wouldn't the operation return zero anyway (unless xsteps==1)? If all we are doing is preventing division by zero errors, then wouldn't we rather want to test for xsteps==1?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with @WeatherGod here, xsteps==1 (and similar for y) is better here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, as a bit of sanity check, is it ever possible for x/ysteps to be zero?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Argument validation for the grid specification arguments could be done on the python side. Thats basically in delaunay.interpolate.LinearInterpolator.getitem (but we should grep for other possible references).
Currently no validation is done. If we want to add them, we should think about the possible semantics of various cases.

To @WeatherGod (2): 0 for x/ysteps would just return an empty grid. If you put a negative value, it would fail before reaching the loop, because numpy raises an exception if you try to allocate an array with negative dimenson.

To @WeatherGod (1), @ianthomas23 : To be exact, I was checking for 0/0 (which is nan) rather than the general */0 (which would be "inf" unless the * is also 0).

The case with xsteps==1 and x0!=x1 may have valid use-cases, but it is not clear where the single point should be. Setting dx=0 in this case is equivalent to putting it on x0. We could, for example place it at 0.5*(x0+x1), or at x1. The stable code would fill the array with inf, which may be considered as some kind of error-indication.
On the other hand, I had actually used the case (x0==x1, xsteps>1) before, and from my point of view I was merely "extending its range of validity" to xsteps==0, without affecting the case described above.

However: I tend to think that having the "*/0" case produce x0 is actually better than the current "inf", and I do agree that it makes the code somewhat more readable.
I'll accept your suggestion - just making sure you understand the implications.

btw: the other fix (#998) has the advantage that the edge cases are naturally resolved, and you do not have to think about all these cases (I guess it would also produce x0 rather than inf, but I did not check that).


rowtri = 0;
for (iy=0; iy<ysteps; iy++) {
Expand Down
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
27 changes: 25 additions & 2 deletions lib/matplotlib/tests/test_delaunay.py
Expand Up @@ -159,7 +159,7 @@ def interpolator(self, func):
z = func(self.x, self.y)
return self.tri.nn_extrapolator(z, bbox=self.xrange+self.yrange)

def make_all_testfuncs(allfuncs=allfuncs):
def make_all_2d_testfuncs(allfuncs=allfuncs):
def make_test(func):
filenames = [
'%s-%s' % (func.func_name, x) for x in
Expand All @@ -186,4 +186,27 @@ def reference_test():
for func in allfuncs:
globals()['test_%s' % func.func_name] = make_test(func)

make_all_testfuncs()
make_all_2d_testfuncs()

# 1d and 0d grid tests

ref_interpolator = Triangulation([0,10,10,0],
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rather than using the module level namespace, it might be worth considering constructing a class which has "test_" methods.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When adding tests to an existing module, I should try staying close to existing style.
make_all_testfuncs() was there (I just renamed it), and since it used module functions rather than class methods, so did I
( and that would have been my default choice anyways in this case).

[0,0,10,10]).linear_interpolator([1,10,5,2.0])

def equal_arrays(a1,a2, tolerance=1e-10):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Were you aware of numpy.testing? In particular the numpy.testing.assert_array_almost_equal function.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nope, but always happy to learn.
help(assert_array_almost_equal) referred me to np.testing.assert_allclose, which in turn referenced np.allclose(). A useful np shortcut I was not aware of. Thanks 👍
(p.s. It seems that generalizations added by assert_allclose, like handling inf/nan values, are redundant in this case, and probably not worth the complication of extra import and runtime, so I use plain allclose() ).

return np.all(np.absolute(a1 - a2) < tolerance)

def test_1d_grid():
res = ref_interpolator[3:6:2j,1:1:1j]
assert equal_arrays(res, [[1.6],[1.9]])

def test_0d_grid():
res = ref_interpolator[3:3:1j,1:1:1j]
assert equal_arrays(res, [[1.6]])

@image_comparison(baseline_images=['delaunay-1d-interp'], extensions=['png'])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should probably remove the fonts from this plot. (see other examples in test_axes)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought the same myself, but wanted to get feedback, since these "freetype_version=" arguments seem to appear everywhere. Indeed, test_axes has a couple of tick-less plots. Point taken.

def test_1d_plots():
x_range = slice(0.25,9.75,20j)
x = np.mgrid[x_range]
for y in xrange(2,10,2):
plt.plot(x, ref_interpolator[x_range,y:y:1j])