[Octave-bug-tracker] [bug #31352] interp2 broken for non-quadratic matrices when using cubic interpolation

[Thorsten Meyer]

Follow-up Comment #15, bug #31352 (project octave):

Hello all,

attached is a changeset containing a merged version of bicubic together with
Marcos changes in interp2 and __splinen__, the tests by Gunnar and a few other
tests. Also, I added to the documentation of the bicubic function.

bicubic will use Marcos algorithm for meshgrid data and my algorithm for
scattered data. Also, there is an undocumented parameter to explicitly choose
which algorithm to use.

As far as I can see, the results are accurate and do not differ by more than
a reasonable numerical roundoff error.

However, as you suspected, Marco, in case of scattered xi and yi, your
algorithm will only be slower than the other one for large x,y or xi, yi
vectors. And I guess it will depend on details of the hardware, where the
breakeven is. I would probably still stick to the simple choice of algorithm
by meshgrid versus scattered. 

Here is a script to compare the performance of the two algorithms for
scattered data and different sizes of x,y,xi,yi:


# compare runtimes of bicubic algorithms
t1 = zeros (30,31);
t2 = zeros (30,31);
o = 1;
more off;
for m = logspace (1, 3.5, 30);
  p = 1;
  for n = logspace (1, 5, 31);
    x  = linspace (-100 * rand(1,1), 100 * rand (1,1), m);
    y  = linspace (-100 * rand(1,1), 100 * rand (1,1), m);
    xi = linspace (min(x), max(x), n);
    yi = linspace (min(x), max(x), n);
    Z = rand (length(y), length(x));

    try
      tic;
      Zi = bicubic(x,y,Z,xi,yi, nan, -0.5, 1);
      t1 (o,p) = toc;
    catch
      t1 (o,p) = nan;
    end_try_catch
#    try
      tic;
      Zi = bicubic(x,y,Z,xi,yi, nan, -0.5, 2);
      t2 (o,p) = toc;
#   catch
#      t2 (o,p) = nan;
#    end_try_catch
 
    p++;
    puts(".");
  endfor
  o++
endfor
more on;

h1=mesh(logspace (1, 5, 31), logspace (1, 3.5, 30), t1);
hold on;
h2=mesh(logspace (1, 5, 31), logspace (1, 3.5, 30), t2);
set(gca, "xscale", "log", "yscale", "log", "zscale", "log")
set(h1, "edgecolor", "blue");
set(h2, "edgecolor", "red");
xlabel ("n xi");
ylabel ("n x");
zlabel ("t(s)");
hold off;


regards

Thorsten



(file #22106)
    _______________________________________________________

Additional Item Attachment:

File name: interp2_fix_caliari            Size:30 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?31352>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/