bug in interpn for method linear and nearest

[Alexander Barth]

On Jan 22, 2008 1:20 AM, David Bateman <adb014 at gmail.com> wrote:
> Alexander Barth wrote:
> >
> > Thank you for reply. I tried to apply your revised patch, but I got this
> > error:
> > laptop:~/Octave/CVS/octave$ patch < interpn.patch
> > patch: **** Overdue `---' at line 16; check line numbers at line 4
> >
> > Do I apply it correctly?
>
> Don't know, and it doesn't matter.. Your patch is better.
>
> > In the attached patch, I added the test-case and removed the commended
> > code as you pointed-out.
> > I also included a test-case for the spline interpolation and we have the
> > same incompatibility:
> > octave:5> [x,y,z] = ndgrid(0:2);
> > octave:6> f = x+y+z;
> > octave:7> fi = interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],'spline')
> > fi =
> >
> > ans(:,:,1) =
> >
> >    1.5000   2.5000
> >    2.5000   3.5000
> >
> > ans(:,:,2) =
> >
> >    2.5000   3.5000
> >    3.5000   4.5000
> >
> > (should be just [1.5000    4.5000])
> >
> > Sorry, I don't understand the spline-code well enough to propose a patch
> > for this.
>
> The spline code treats the interpolation in each dimension separately
> with the 1D spline code, It iteratively generates the yi value over the
> dimensions, and I don't see an easy way of getting rid of the ndgridding
> as the ndgridded data is needed at the intermediate stages of the
> calculation. The only way I see to get rid of it is with a patch like
> the attached, that implements a sort of ND diag function.
>
> However, I wonder if the patch you propose is completely correct.
> Consider the case
>
> [x,y,z] = ndgrid(0:2);
> f = x+y+z;
> interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5 2])
>
> The size of the yi variables are not equal and I suspect (though can't
> test for lack of matlab at the moment) that in that case it should be
> ndgridded. The attached patch also implements that behavior, though this
> should be checked before the patch is applied.
>
> D.
>

Matlab returns in this case an error:
>> [x,y,z] = ndgrid(0:2);
>> f = x+y+z;
>> interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5 2])
??? Error using ==> interpn
Y1,Y2,Y3, etc. must all the same size or vectors of different orientations.

The patched octave version too:
octave:6> interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5 2])
error: interpn: incompatible size of argument number 7

As I interprete matlab's error message (vectors of different
orientations), I expected the following to be ndgridded:

>> interpn(x,y,z,f,reshape([.5 1.5],[2 1 1]),reshape([.5 1.5],[1 2
1]),reshape([.5 1.5 2],[1 1 3]))
??? Error using ==> interpn
Y1,Y2,Y3, etc. must all the same size or vectors of different orientations.

But this seems not to be the case. As far as I have tested, we are
compatible relative to the acceptable input parameters (and we don't
produce a confusing error messages!).

>
> As the underlying spline function has this behavior of ndgridding the
> output values if yi is a vector and the values array is an NDArray and
> does in both Octave and Matlab, I'm not sure the right way to fix this.
> A quick hack would be to do the equivalent of an ND diag function to get
> the required values.
>

I agree with that approach; to make the spline interpolation
compatible now and to worry about the optimization later. Do we
already have already an equivalent of an ND diag function?

Alex