bugGNU Octave - Bugs: bug #45542, griddatan unable to perform...

 
 

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

Add a New Comment Rich Markup
   

Jump to the original submission

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>
Group administrator
Fri 12 Jul 2019 05:26:43 PM UTC, comment #7: 

I'll take a look.

Rik <rik5>
Group administrator
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>
Group administrator
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>
Group administrator
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>
Group administrator
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>
Group Member
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.

(file #34467)

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

 

(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)

 

Depends on the following items: None found

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 group members can vote.

     

    Follow 17 latest changes.

    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
        Attached File- Added test_griddatan.mat, #34437

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code