speed of interp2

[Peter Cloetens]

Hello,

It seems that the function interp2 got much slower at a certain point
(~factor 6 for linear interpolation and ~factor 40 for nearest neighbour).
We use octave version 2.1.73.
The old interp2 is older than 2005-04-23 but prior to 2005-09-07.
The new interp2 corresponds to octave-forge-2006.03.17 and it is also the one included in recent versions of octave.

Reading the history and comparing the dates of the two versions of
interp2 I suspect that it is linked to the modification:
(XI,YI no longer need to be "meshgridded")

Here are some examples:
old interp2:
nearest:
> z=rand(1024);tic;zi=interp2_old(z,1:0.5:1024,(1:0.5:1024)','nearest');toc
ans = 0.43838
linear:
> z=rand(1024);tic;zi=interp2_old(z,1:0.5:1024,(1:0.5:1024)','linear');toc
ans = 4.1426
no bicubic available

recent interp2:
nearest:
> z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','nearest');toc
ans = 18.197
linear:
> z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','linear');toc
ans = 26.559
cubic:
> z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','cubic');toc
ans = 33.532
cubic calling bicubic directly:
> z=rand(1024);tic;zi=bicubic(1:1024,(1:1024)',z,1:0.5:1024,(1:0.5:1024)');toc
ans = 10.435

In the old code one has:
    xidx = lookup(X(1, 2:end-1), XI(1,:))' + 1;
so xidx will have length 2047

In the recent code one has
    xidx = lookup(X(1, 2:end-1), XI) + 1;
so xidx will have length 2047*2047

Maybe a separate code should be foreseen in the case XI and YI are vectors.

Cheers,
Peter


-------- Forwarded Message --------

From: Peter Cloetens <Peter%20Cloetens%20%3ccloetens@...>
To: cloetens@...
Date: Sun, 4 Feb 2007 10:43:42 +0100

Begin forwarded message:

From: David Bateman <David.Bateman@...>
Date: January 31, 2007 4:00:27 PM GMT+01:00
To: Peter Cloetens <cloetens@...>
Cc: pkienzle@..., balden@...
Subject: Re: speed of interp2

Peter Cloetens wrote:
Hello,
We are updating our octave-forge version.
It seems that the function interp2 got much slower at a certain point
(~factor 5).
Reading the history and comparing the dates of the two versions of
interp2 I suspect that it is linked to the modification:
(XI,YI no longer need to be "meshgridded")
Fortunately it is possible to use the old version of interp2 ;-)
Peter


## Author: Kai Habel <kai.habel@... <mailto:kai.habel@...>>
## 2005-03-02 Thomas Weber <weber@...
<mailto:weber@...>>
##     * Add test cases
## 2005-03-02 Paul Kienzle <pkienzle@...
<mailto:pkienzle@...>>
##     * Simplify
## 2005-04-23 Dmitri A. Sergatskov <dasergatskov@...
<mailto:dasergatskov@...>>
##     * Modified demo and test for new gnuplot interface
## 2005-09-07 Hoxide <hoxide_dirac@...
<mailto:hoxide_dirac@...>>
##     * Add bicubic interpolation method
##     * Fix the eat line bug when the last element of XI or YI is
negative or zero.
## 2005-11-26 Pierre Baldensperger <balden@...
<mailto:balden@...>>
##     * Rather big modification (XI,YI no longer need to be
##       "meshgridded") to be consistent with the help message
##       above and for compatibility.

Peter,

This not my code and so it will be difficult for me to debug the issue
you raise. However, if you supply a small benchmark test of interp2 and
the CVS dates of the versions of interp2 that seemed to work and those
that didn't, I can attempt to find the issue...

I took a quick look, and its seems in the simple tests I tried that it
wasn't the gridding that took the largest portion of the time. Also
which type of interpolation is slow? Is it the "linear" interpolation
that is the default, one of the others, or all of them?

Note also that interp2 is now part of Octave itself and so this might be
better to be discussed on bug@......

Regards
David

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On  4-Feb-2007, Peter Cloetens wrote:

| Hello,
|
| It seems that the function interp2 got much slower at a certain point
| (~factor 6 for linear interpolation and ~factor 40 for nearest
| neighbour).
| We use octave version 2.1.73.
| The old interp2 is older than 2005-04-23 but prior to 2005-09-07.
| The new interp2 corresponds to octave-forge-2006.03.17 and it is also
| the one included in recent versions of octave.
|
| Reading the history and comparing the dates of the two versions of
| interp2 I suspect that it is linked to the modification:
| (XI,YI no longer need to be "meshgridded")
|
| Here are some examples:
| old interp2:
| nearest:
| > z=rand(1024);tic;zi=interp2_old(z,1:0.5:1024,
| (1:0.5:1024)','nearest');toc
| ans = 0.43838
| linear:
| > z=rand(1024);tic;zi=interp2_old(z,1:0.5:1024,
| (1:0.5:1024)','linear');toc
| ans = 4.1426
| no bicubic available
|
| recent interp2:
| nearest:
| > z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','nearest');toc
| ans = 18.197
| linear:
| > z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','linear');toc
| ans = 26.559
| cubic:
| > z=rand(1024);tic;zi=interp2(z,1:0.5:1024,(1:0.5:1024)','cubic');toc
| ans = 33.532
| cubic calling bicubic directly:
| > z=rand(1024);tic;zi=bicubic(1:1024,(1:1024)',z,1:0.5:1024,
| (1:0.5:1024)');toc
| ans = 10.435
|
| In the old code one has:
|     xidx = lookup(X(1, 2:end-1), XI(1,:))' + 1;
| so xidx will have length 2047
|
| In the recent code one has
|     xidx = lookup(X(1, 2:end-1), XI) + 1;
| so xidx will have length 2047*2047
|  
| Maybe a separate code should be foreseen in the case XI and YI are
| vectors.

If you have a fix, how about submitting a patch?

Thanks,

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Peter Cloetens]

John W. Eaton wrote:

On  4-Feb-2007, Peter Cloetens wrote:

| Hello,
|
| It seems that the function interp2 got much slower at a certain point
| (~factor 6 for linear interpolation and ~factor 40 for nearest
| neighbour).
|...
| Reading the history and comparing the dates of the two versions of
| interp2 I suspect that it is linked to the modification:
| (XI,YI no longer need to be "meshgridded")
|...
| Maybe a separate code should be foreseen in the case XI and YI are
| vectors.

If you have a fix, how about submitting a patch?

Thanks,

jwe
_______________________________________________
Bug-octave mailing list
[e-mail]
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

Hello,
I would like to come back to the issue of the slow interp2. As it is used by other functions (imrotate, imremap, ...) this is quite important.
I made a faster version:
interp2.m

It uses faster code when XI and YI are vectors and when X and Y are not given. Also in the general case it is faster.
It also removes a bug in the cubic interpolation that would give an error if the matrix is not square.

I made a simple test program:
interp2_benchmark.m

The result with the old code:
cloetens@coral39:cloetens> interp2_benchmark(2048)
Interpolating 2047 x 2048 matrix ; nearest ; XI and YI are matrix : 11.6679 seconds
Interpolating 2047 x 2048 matrix ; nearest ; XI and YI are vector : 11.3098 seconds
Interpolating 2047 x 2048 matrix ; nearest ; X, Y, XI and YI are matrix : 11.1117 seconds
Interpolating 2047 x 2048 matrix ; linear ; XI and YI are matrix : 14.8497 seconds
Interpolating 2047 x 2048 matrix ; linear ; XI and YI are vector : 14.8174 seconds
Interpolating 2047 x 2048 matrix ; linear ; X, Y, XI and YI are matrix : 14.6365 seconds
error: X and Y size must match Z dimensions
error: evaluating if command near line 246, column 7
error: evaluating if command near line 243, column 5
error: evaluating if command near line 164, column 3
error: called from `interp2' in file `/sware/pub/octave-3.0.1/redhate4.x86_64/share/octave/3.0.1/m/general/interp2.m'
error: evaluating for command near line 33, column 5
error: called from `interp2_benchmark' in file `/users/cloetens/octave/interp2_benchmark.m'

The result with the new code:
cloetens@coral39:cloetens> interp2_benchmark(2048)
Interpolating 2047 x 2048 matrix ; nearest ; XI and YI are matrix : 2.31001 seconds
Interpolating 2047 x 2048 matrix ; nearest ; XI and YI are vector : 0.333178 seconds
Interpolating 2047 x 2048 matrix ; nearest ; X, Y, XI and YI are matrix : 7.03896 seconds
Interpolating 2047 x 2048 matrix ; linear ; XI and YI are matrix : 5.70732 seconds
Interpolating 2047 x 2048 matrix ; linear ; XI and YI are vector : 2.33167 seconds
Interpolating 2047 x 2048 matrix ; linear ; X, Y, XI and YI are matrix : 10.4666 seconds
Interpolating 2047 x 2048 matrix ; cubic ; XI and YI are matrix : 5.01822 seconds
Interpolating 2047 x 2048 matrix ; cubic ; XI and YI are vector : 4.8136 seconds
Interpolating 2047 x 2048 matrix ; cubic ; X, Y, XI and YI are matrix : 4.38982 seconds

Peter

[John W. Eaton]

On 26-Oct-2008, Peter Cloetens wrote:

|
|
| John W. Eaton wrote:
| >
| > On  4-Feb-2007, Peter Cloetens wrote:
| >
| > | Hello,
| > |
| > | It seems that the function interp2 got much slower at a certain point
| > | (~factor 6 for linear interpolation and ~factor 40 for nearest
| > | neighbour).
| > |...
| > | Reading the history and comparing the dates of the two versions of
| > | interp2 I suspect that it is linked to the modification:
| > | (XI,YI no longer need to be "meshgridded")
| > |...
| > | Maybe a separate code should be foreseen in the case XI and YI are
| > | vectors.
| >
| > If you have a fix, how about submitting a patch?
| >
| > Thanks,
| >
| > jwe
| > _______________________________________________
| > Bug-octave mailing list
| > Bug-octave@...
| > https://www.cae.wisc.edu/mailman/listinfo/bug-octave
| >
| >
|
| Hello,
| I would like to come back to the issue of the slow interp2. As it is used by
| other functions (imrotate, imremap, ...) this is quite important.
| I made a faster version:
| http://www.nabble.com/file/p20177305/interp2.m interp2.m

Will you please submit a patch that is relative to the version you
started with?  That way we can see precisely what changes you made,
and we will also avoid possibly wiping out other changes that might
have been made in Octave's sources in the time since the version you
started with was released.

jw
e
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Peter Cloetens]

John W. Eaton wrote:

> On 26-Oct-2008, Peter Cloetens wrote:
>
> |
> |
> | John W. Eaton wrote:
> | >
> | > On  4-Feb-2007, Peter Cloetens wrote:
> | >
> | > | Hello,
> | > |
> | > | It seems that the function interp2 got much slower at a certain point
> | > | (~factor 6 for linear interpolation and ~factor 40 for nearest
> | > | neighbour).
> | > |...
> | > | Reading the history and comparing the dates of the two versions of
> | > | interp2 I suspect that it is linked to the modification:
> | > | (XI,YI no longer need to be "meshgridded")
> | > |...
> | > | Maybe a separate code should be foreseen in the case XI and YI are
> | > | vectors.
> | >
> | > If you have a fix, how about submitting a patch?
> | >
> | > Thanks,
> | >
> | > jwe
> | > _______________________________________________
> | > Bug-octave mailing list
> | > Bug-octave@...
> | > https://www.cae.wisc.edu/mailman/listinfo/bug-octave
> | >
> | >
> |
> | Hello,
> | I would like to come back to the issue of the slow interp2. As it is used by
> | other functions (imrotate, imremap, ...) this is quite important.
> | I made a faster version:
> | http://www.nabble.com/file/p20177305/interp2.m interp2.m
>
> Will you please submit a patch that is relative to the version you
> started with?  That way we can see precisely what changes you made,
> and we will also avoid possibly wiping out other changes that might
> have been made in Octave's sources in the time since the version you
> started with was released.
>
> jw
> e

I'm not familiar with creating a patch, but I will try to find out how
to do this. I started with the version of interp2.m that comes with
Octave 3.0.3, but the optimized code for XI and YI vectors is indeed
older. The versions passes all tests and demos successfully, but I will
check if there were intermediate changes that I would have missed.
Also, I think more test cases should be added especially for the cubic
interpolation.
Peter
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Søren Hauberg]

tir, 28 10 2008 kl. 14:47 +0100, skrev Peter Cloetens:
> I'm not familiar with creating a patch, but I will try to find out how
> to do this. I started with the version of interp2.m that comes with
> Octave 3.0.3,

If you use Linux or Unix, creating a patch is fairly easy. In a terminal
simple type

  diff -u interp2.m your_new_interp2.m > interp2.patch

here 'interp2.m' should point to the version from 3.0.3,
'your_new_interp2.m' should be your file, and the resulting patch will
be stored in the file 'interp2.patch'. Of course, it's better if you use
Mercurial to do this, but the way I've just described should be just
fine.

Søren

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Sergei Steshenko-2]

--- On Tue, 10/28/08, Peter Cloetens <cloetens@...> wrote:

> From: Peter Cloetens <cloetens@...>
> Subject: Re: speed of interp2 and bug with method cubic
> To: bug-octave@...
> Date: Tuesday, October 28, 2008, 6:47 AM
> John W. Eaton wrote:
> > On 26-Oct-2008, Peter Cloetens wrote:
> >
> > |
> > |
> > | John W. Eaton wrote:
> > | >
> > | > On  4-Feb-2007, Peter Cloetens wrote:
> > | >
> > | > | Hello,
> > | > |
> > | > | It seems that the function interp2 got much
> slower at a certain point
> > | > | (~factor 6 for linear interpolation and
> ~factor 40 for nearest
> > | > | neighbour).
> > | > |...
> > | > | Reading the history and comparing the dates
> of the two versions of
> > | > | interp2 I suspect that it is linked to the
> modification:
> > | > | (XI,YI no longer need to be
> "meshgridded")
> > | > |...
> > | > | Maybe a separate code should be foreseen in
> the case XI and YI are
> > | > | vectors.
> > | >
> > | > If you have a fix, how about submitting a
> patch?
> > | >
> > | > Thanks,
> > | >
> > | > jwe
> > | > _______________________________________________
> > | > Bug-octave mailing list
> > | > Bug-octave@...
> > | >
> https://www.cae.wisc.edu/mailman/listinfo/bug-octave
> > | >
> > | >
> > |
> > | Hello,
> > | I would like to come back to the issue of the slow
> interp2. As it is used by
> > | other functions (imrotate, imremap, ...) this is
> quite important.
> > | I made a faster version:
> > | http://www.nabble.com/file/p20177305/interp2.m
> interp2.m
> >
> > Will you please submit a patch that is relative to the
> version you
> > started with?  That way we can see precisely what
> changes you made,
> > and we will also avoid possibly wiping out other
> changes that might
> > have been made in Octave's sources in the time
> since the version you
> > started with was released.
> >
> > jw
> > e
> I'm not familiar with creating a patch, but I will try
> to find out how
> to do this. I started with the version of interp2.m that
> comes with
> Octave 3.0.3, but the optimized code for XI and YI vectors
> is indeed
> older. The versions passes all tests and demos
> successfully, but I will
> check if there were intermediate changes that I would have
> missed.
> Also, I think more test cases should be added especially
> for the cubic
> interpolation.
> Peter
> _______________________________________________
> Bug-octave mailing list
> Bug-octave@...
> https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

Peter, 'man patch' (in the end) gives an easy to remember example:

"
NOTES FOR PATCH SENDERS
       There are several things you should bear in mind if you are going to be sending out patches.

       Create  your  patch  systematically.  A good method is the command diff -Naur old new where old and new identify the old and new directo‐
       ries.  The names old and new should not contain any slashes.  The diff command's headers should have dates and times  in  Universal  Time
       using  traditional  Unix  format,  so that patch recipients can use the -Z or --set-utc option.  Here is an example command, using Bourne
       shell syntax:

          LC_ALL=C TZ=UTC0 diff -Naur gcc-2.7 gcc-2.8
".

However, taking into account the possibility of existence of outdated code
in the file with your changes you probably better send the file as is.

If you are not familiar with 'tkdiff' - try it - makes it easier to
comprehend differences. There are also some 3 way diff utilities, like
KDiff3.


Regards,
  Sergei.


     

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton-3]

On 26-Oct-2008, Peter Cloetens wrote:

| John W. Eaton wrote:
| >
| > On  4-Feb-2007, Peter Cloetens wrote:
| >
| > | Hello,
| > |
| > | It seems that the function interp2 got much slower at a certain point
| > | (~factor 6 for linear interpolation and ~factor 40 for nearest
| > | neighbour).
| > |...
| > | Reading the history and comparing the dates of the two versions of
| > | interp2 I suspect that it is linked to the modification:
| > | (XI,YI no longer need to be "meshgridded")
| > |...
| > | Maybe a separate code should be foreseen in the case XI and YI are
| > | vectors.
| >
| > If you have a fix, how about submitting a patch?
| >
| > Thanks,
| >
| > jwe
| > _______________________________________________
| > Bug-octave mailing list
| > Bug-octave@...
| > https://www.cae.wisc.edu/mailman/listinfo/bug-octave
| >
| >
|
| Hello,
| I would like to come back to the issue of the slow interp2. As it is used by
| other functions (imrotate, imremap, ...) this is quite important.
| I made a faster version:
| http://www.nabble.com/file/p20177305/interp2.m interp2.m
|
| It uses faster code when XI and YI are vectors and when X and Y are not
| given. Also in the general case it is faster.
| It also removes a bug in the cubic interpolation that would give an error if
| the matrix is not square.

It doesn't look like these changes were ever incorporated.  Could
someone who is more familiar with the interpolation functions take a
look and comment?  Also, it would be helpful if someone could generate
a patch.  Peter, what vesion of Octave was your modified interp2 based
on?  Having a diff relative to that would be useful.

| ## 2008-10-25 Peter Cloetens <cloetens@...>
| ##     * use old, faster code when XI and YI are given as vectors
| ##     * make code faster ; optimized when X, Y are default
| ##     * fix problem in cubic interpolation when Z is not square

Please don't add this information to the files.  ChangeLog info
belongs in the ChangeLog files.  I know that there is info like this
in some files in Octave, usually because these lines were left over
from Octave Forge.  What should we do with them?  I don't want to
encourage cluttering all the files in Octave with ChangeLog info.

Thanks,

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave