akima 0.4-5, interpp() bug = COMMON block problem

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

akima 0.4-5, interpp() bug = COMMON block problem

Albrecht Gebhardt
Hi,

I'm currently hunting a bug in the akima library, especially in the code
behind the interpp.old function (bi-variate linear interpolation).
It is based on a triangulation algorithm, interpolation at a given point
needs to know the triangle which contains this point, then the
interpolation is a straightforward calculation based on the three
vertexes.

The problem is: Sometimes the triangle indices are not given back
correctly, they just default to 1 leading to wrong results.

The following lines can be used to visualize it using the rgl library.
If the error occurs (may be architecture depending) the interpolation
"steepest" part of the interp() surface will not be hit by the interpp()
interpolation points.

library(akima)
library(rgl)
data(akima)
# data
rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red")
rgl.bbox()
# interp:
akima.li <- interp(akima$x, akima$y, akima$z,
                   xo=seq(min(akima$x), max(akima$x), length = 200),
                   yo=seq(min(akima$y), max(akima$y), length = 200))
# interp surface:
rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
# interpp:
akima.p <- interpp(akima$x, akima$y, akima$z,
                    runif(2000,0,25),
                    runif(2000,0,20))
# interpp (partially wrong) points:
rgl.points(akima.p$x,akima.p$z , akima.p$y,size=4,color="yellow")

The errors occurs at least with R 2.1.1 on a i386 Ubuntu Breezy system.
If I compile without "-O2" in FFLAGS the error vanishes.

Meanwhile I could track down the bug to the use of COMMON blocks in the
underlying Fortran code:

During the interpp() call a Fortran routine (IDLCTN called from IDBVIP)
is called several times. Only during the first call a COMMON block is
initialized, subsequent calls rely on the values initialized by the
first call. But these values just vanish after the first call.

I already have workaround: I call the initialization on every entry,
this fixes the problem.

===================================================================
RCS file: /home/cvs/math/stat/R/akima/src/idlctn.f,v
retrieving revision 1.2
diff -u -r1.2 idlctn.f
--- idlctn.f    19 Aug 1998 21:54:17 -0000      1.2
+++ idlctn.f    1 Feb 2006 10:28:53 -0000
@@ -49,7 +49,7 @@
       X0 = XII
       Y0 = YII
 C PROCESSING FOR A NEW SET OF DATA POINTS
-      IF (NIT.NE.0) GO TO 80
+C      IF (NIT.NE.0) GO TO 80
       NIT = 1
 C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS.
       XMN = XD(1)


I guess, using COMMON blocks is generally a bad idea when using Fortran
code within R? Shall I leave it with my workaround or should I search
for more details of the COMMON block misuse?



Best regards

Albrecht Gebhardt


--
// Albrecht Gebhardt          Tel.: (++43 463) 2700/3118
// Institut fuer Mathematik   Fax : (++43 463) 2700/3198
// Universitaet Klagenfurt    mailto:[hidden email]
// Universitaetsstr. 65       http://www.math.uni-klu.ac.at/~agebhard
// A-9020 Klagenfurt, Austria
// GPG PK: http://www.math.uni-klu.ac.at/~agebhard/agebhard.asc
// GPG FP: F46F 656E E83C 9323 CE30  FF8F 9DBA D1A3 B55A 78A6

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: akima 0.4-5, interpp() bug = COMMON block problem

Albrecht Gebhardt
After reading my own question once again, I think I can answer it
myself:

The code had too _few_ COMMON blocks. It seems the code relied on that,
that several local variables would carry their values after being
initialised in the first subroutine call until the second call. This is
not the case. I introduced a new COMMON block containing (at least) all
these variables and this fixes it.

After I have done some more checks, you should expect a new version of
akima soon.


Best wishes

Albrecht


Am Mittwoch, den 01.02.2006, 11:31 +0100 schrieb Albrecht Gebhardt:

> Hi,
>
> I'm currently hunting a bug in the akima package, especially in the code
> behind the interpp.old function (bi-variate linear interpolation).
> It is based on a triangulation algorithm, interpolation at a given point
> needs to know the triangle which contains this point, then the
> interpolation is a straightforward calculation based on the three
> vertexes.
>
> The problem is: Sometimes the triangle indices are not given back
> correctly, they just default to 1 leading to wrong results.
>
> The following lines can be used to visualize it using the rgl package.
> If the error occurs (may be architecture depending) the interpolation
> "steepest" part of the interp() surface will not be hit by the interpp()
> interpolation points.
>
> library(akima)
> library(rgl)
> data(akima)
> # data
> rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red")
> rgl.bbox()
> # interp:
> akima.li <- interp(akima$x, akima$y, akima$z,
>                    xo=seq(min(akima$x), max(akima$x), length = 200),
>                    yo=seq(min(akima$y), max(akima$y), length = 200))
> # interp surface:
> rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
> # interpp:
> akima.p <- interpp(akima$x, akima$y, akima$z,
>                     runif(2000,0,25),
>                     runif(2000,0,20))
> # interpp (partially wrong) points:
> rgl.points(akima.p$x,akima.p$z , akima.p$y,size=4,color="yellow")
>
> The errors occurs at least with R 2.1.1 on a i386 Ubuntu Breezy system.
> If I compile without "-O2" in FFLAGS the error vanishes.
>
> Meanwhile I could track down the bug to the use of COMMON blocks in the
> underlying Fortran code:
>
> During the interpp() call a Fortran routine (IDLCTN called from IDBVIP)
> is called several times. Only during the first call a COMMON block is
> initialized, subsequent calls rely on the values initialized by the
> first call. But these values just vanish after the first call.
>
> I already have workaround: I call the initialization on every entry,
> this fixes the problem.
>
> ===================================================================
> RCS file: /home/cvs/math/stat/R/akima/src/idlctn.f,v
> retrieving revision 1.2
> diff -u -r1.2 idlctn.f
> --- idlctn.f    19 Aug 1998 21:54:17 -0000      1.2
> +++ idlctn.f    1 Feb 2006 10:28:53 -0000
> @@ -49,7 +49,7 @@
>        X0 = XII
>        Y0 = YII
>  C PROCESSING FOR A NEW SET OF DATA POINTS
> -      IF (NIT.NE.0) GO TO 80
> +C      IF (NIT.NE.0) GO TO 80
>        NIT = 1
>  C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS.
>        XMN = XD(1)
>
>
> I guess, using COMMON blocks is generally a bad idea when using Fortran
> code within R? Shall I leave it with my workaround or should I search
> for more details of the COMMON block misuse?
>
>
>
> Best regards
>
> Albrecht Gebhardt
>
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
--
// Albrecht Gebhardt          Tel.: (++43 463) 2700/3118
// Institut fuer Mathematik   Fax : (++43 463) 2700/3198
// Universitaet Klagenfurt    mailto:[hidden email]
// Universitaetsstr. 65       http://www.math.uni-klu.ac.at/~agebhard
// A-9020 Klagenfurt, Austria
// GPG PK: http://www.math.uni-klu.ac.at/~agebhard/agebhard.asc
// GPG FP: F46F 656E E83C 9323 CE30  FF8F 9DBA D1A3 B55A 78A6

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

signature.asc (196 bytes) Download Attachment