# bug #45542: griddatan unable to perform 4-dimensional interpolation.

 Submitter: None Submitted: Mon 13 Jul 2015 02:25:03 PM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Incorrect Result Status: Fixed Assigned to: None Originator Name: Juho Iipponen Originator Email: -email is unavailable- Open/Closed:  * Closed Release:  * dev Operating System:  * Any Fixed Release: None Planned Release: None
* Mandatory Fields

Comment Type & Canned Response:

## Discussion

Fri 12 Jul 2019 07:43:04 PM UTC, comment #8:

I checked in a new fix here (https://hg.savannah.gnu.org/hgweb/octave/rev/3f354ef16400).

Apparently the code had not been reviewed in a very long while because there were lots of obvious problems.  We were calculated an intermediate variable three times rather than just once and re-using it.  And there was this bit of code that was not used at all, but involved quite a bit of calculation.

-    ## assign x,y for each point of simplex
-    xt = reshape (x(tri(tri_list,:),:), [nr_t, n+1, n]);
-    yt = y(tri(tri_list,:));

The problem turned out to be caused by the orientation rules when doing vector/matrix indexing versus vector/vector indexing.  In the first, the result takes the shape of the matrix index.  In the latter, the shape is that of the vector being indexed rather than the index vector.  I added a special case to resolve this.

+    if (isscalar (tri_list))
+      ## Special case required by orientation rules for vector/vector index.
+      yi(valid) = sum (y(tri(tri_list,:)).' .* bary_list, 2);
+    else
+      yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2);
endif

Rik <rik5>
Fri 12 Jul 2019 05:26:43 PM UTC, comment #7:

I'll take a look.

Rik <rik5>
Fri 12 Jul 2019 04:10:59 PM UTC, comment #6:

Yes, I'm seeing the same three test failures after changeset bee80e27dcc5 and none before.  I haven't checked to see, but I'm guessing that the tests need to be adjusted to match the new expected behavior?

John W. Eaton <jwe>
Fri 12 Jul 2019 11:59:14 AM UTC, comment #5:

The test suite running on the build bots fails with the following errors:

>>>>> processing /scratch/buildbot/workers/jwe-debian-x86_64-1/clang-4_0-debian/src/scripts/geometry/griddata3.m
***** testif HAVE_QHULL
old_state = rand ("state");
restore_state = onCleanup (@() rand ("state", old_state));
rand ("state", 0);
x = 2 * rand (1000, 1) - 1;
y = 2 * rand (1000, 1) - 1;
z = 2 * rand (1000, 1) - 1;
v = x.^2 + y.^2 + z.^2;
[xi, yi, zi] = meshgrid (-0.8:0.2:0.8);
vi = griddata3 (x, y, z, v, xi, yi, zi, "linear");
vv = vi - xi.^2 - yi.^2 - zi.^2;
assert (max (abs (vv(:))), 0, 0.1);
!!!!! test failed
product: nonconformant arguments (op1 is 4x729, op2 is 729x4)
***** testif HAVE_QHULL
old_state = rand ("state");
restore_state = onCleanup (@() rand ("state", old_state));
rand ("state", 0);
x = 2 * rand (1000, 1) - 1;
y = 2 * rand (1000, 1) - 1;
z = 2 * rand (1000, 1) - 1;
v = x.^2 + y.^2 + z.^2;
[xi, yi, zi] = meshgrid (-0.8:0.2:0.8);
vi = griddata3 (x, y, z, v, xi, yi, zi, "nearest");
vv = vi - xi.^2 - yi.^2 - zi.^2;
assert (max (abs (vv(:))), 0, 0.1);
!!!!! test failed
product: nonconformant arguments (op1 is 4x729, op2 is 729x4)
>>>>> processing /scratch/buildbot/workers/jwe-debian-x86_64-1/clang-4_0-debian/src/scripts/geometry/griddatan.m
***** testif HAVE_QHULL
[xx,yy] = meshgrid (linspace (-1,1,32));
xi = [xx(:), yy(:)];
x = 2*rand (100,2) - 1;
x = [x;1,1;1,-1;-1,-1;-1,1];
y = sin (2 * sum (x.^2,2));
zz = griddatan (x,y,xi,"linear");
zz2 = griddata (x(:,1),x(:,2),y,xi(:,1),xi(:,2),"linear");
assert (zz, zz2, 1e-10);
!!!!! test failed
product: nonconformant arguments (op1 is 3x1024, op2 is 1024x3)

Re-opening because I think this might be related to this change.

Markus Mützel <mmuetzel>
Fri 12 Jul 2019 03:18:28 AM UTC, comment #4:

I checked in the change suggested by Kai under Juho's name here (https://hg.savannah.gnu.org/hgweb/octave/rev/bee80e27dcc5).

Marking as fixed and closing report.

Rik <rik5>
Tue 18 Jun 2019 07:17:41 AM UTC, comment #3:

The patch of comment #1 (only three lines changed in file #34467) proved as useful for bug #56515.  Here a proper diff with Octave code style improvements:

diff -r f7bf27159cf7 scripts/geometry/griddatan.m
--- a/scripts/geometry/griddatan.m        Fri Jun 14 10:32:06 2019 -0500
+++ b/scripts/geometry/griddatan.m        Tue Jun 18 09:10:34 2019 +0200
@@ -81,7 +81,9 @@ function yi = griddatan (x, y, xi, metho
yt = y(tri(tri_list,:));

## Use barycentric coordinate of point to calculate yi
-    yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2);
+    if (! isempty (bary_list))
+      yi(valid) = sum (y(tri(tri_list,:))' .* bary_list, 2);
+    endif

else
error ("griddatan: unknown interpolation METHOD");

I can confirm, that this issue is still present in Octave 5.1.0 and works with Matlab R2019a.

Thus I suggest to merge this patch to stable with a proper test credited to you "Juho Iipponen <juho.iipponen@helsinki.fi>", if you agree.

Kai Torben Ohlhus <siko1056>
Sat 02 Jul 2016 02:08:24 AM UTC, comment #2:

Thanks for the patch, Juho.

Could you please prepare a diff of the file from the original, rather than providing a complete new file?  That makes it easier to apply to the current codebase.

Lachlan Andrew <lachlan>
Mon 20 Jul 2015 01:11:15 PM UTC, comment #1:

Taking the transpose of "y(tri(tri_list,:))" before multiplication (see original submission) seems to be the correct course of action. However, it only works if a point lies within the convex hull, when both vectors to be multiplied are non-empty. Therefore, I have added a guard before that line:

if ~isempty(bary_list)
yi(valid) = sum (y(tri(tri_list,:))' .* bary_list, 2);
end

I have attached the fixed griddatan.m function, but I would want someone more knowledgeable to have a look at it.

Juho Iipponen <juppiega>
Mon 13 Jul 2015 02:25:03 PM UTC, original submission:

The griddatan function seems to be unable to perform 4-dimensional interpolation (in e.g. (x, y, z, time)-coordinates). While running the attached script "test_griddatan.m" with test input data "test_griddatan.mat", the function fails while attempting to fill in the interpolated values for output (last line of griddatan.m to be executed for linear interpolation). The error message is:

error: griddatan: A(I) = X: X must have the same size as I
error: called from
griddatan at line 84 column 15
test_griddatan at line 2 column 26

To me, this seems to be caused by different shape of vectors "y(tri(tri_list,:))" and "bary_list", which are multiplied element-wise. In this example, the former is of size 5x1, while the latter is of size 1x5. If I transpose "y(tri(tri_list,:))" before multiplication, the last line yields a result of 3.1866e-14, which is very close to the one given by Matlab (3.1858e-14). However, I have no experience with multidimensional geometry, so I'm not sure if taking the transpose is the correct solution.

Anonymous

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #34467:  griddatan.m added by juppiega (4KiB - text/x-objcsrc)
file #34436:  test_griddatan.m added by None (106B - text/x-objcsrc)
file #34437:  test_griddatan.mat added by None (1KiB - application/octet-stream)

## Dependencies

Depends on the following items: None found

Items that depend on this one

Digest:
bug dependencies.

Carbon-Copy List
• -email is unavailable- added by jwe (Posted a comment)
• -email is unavailable- added by mmuetzel (Posted a comment)
• -email is unavailable- added by siko1056 (Posted a comment)
• -email is unavailable- added by lachlan (Posted a comment)
• -email is unavailable- added by rik5 (Updated the item)
• -email is unavailable- added by juppiega (Updated the item)
• -email is unavailable- added by None (Submitted the item)
•

There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

Only project members can vote.

## History

Date Changed by Updated Field Previous Value => Replaced by
2019-07-12 rik5 Open/ClosedOpen Closed
2019-07-12 rik5 StatusIn Progress Fixed
2019-07-12 mmuetzel StatusFixed In Progress
Open/ClosedClosed Open
2019-07-12 rik5 StatusPatch Submitted Fixed
Open/ClosedOpen Closed
Release5.1.0 dev
2019-06-18 siko1056 StatusIn Progress Patch Submitted
Release4.0.0 5.1.0
Operating SystemGNU/Linux Any
2019-06-18 siko1056 Dependencies- bugs #56515 is dependent
2016-07-02 lachlan StatusPatch Submitted In Progress
2015-12-01 rik5 Item GroupSegfault, Bus Error, etc. Incorrect Result
StatusNone Patch Submitted
2015-07-20 juppiega Attached File- Added griddatan.m, #34467
2015-07-13 None Attached File- Added test_griddatan.m, #34436