R devel - Jan 2012 - Numerical instability in new R Windows development version

I have a question concerning the new Windows toolchain for R >= 2.14.2.
When trying out my package 'pracma' on the win-builder development
version
it will stop with the following error message:

  > f3 <- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1))
  > dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err = 2e-8
  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)) : NaNs produced
  Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)) : NaNs produced
  Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs,  :
    non-finite function value
  Calls: dblquad ...
         <Anonymous> -> f -> do.call -> mapply ->
<Anonymous> -> integrate
  Execution halted
  ** running examples for arch 'x64' ... ERROR
  Running examples in 'pracma-Ex.R' failed

This probably means that the following expression got negative for some
values x, y:

  (1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)

It appears to be an often used trick in numerical analysis. One advantage is
that a function using it is immediately vectorized while an expression such
as, e.g., "max(0, 1 - (x^2 + y^2))" is not.

The example runs fine on Debian Linux and Mac OS X 32-/64-bit architectures.
In my understanding the approach is correct and, as said above, often used in
numerical applications.

Can someone explain to me why this fails for the Windows 64-bit compiler and
what I should use instead. Thanks.

Hans Werner Borchers
ABB Corporate Research

On 12-01-27 7:23 AM, Hans W Borchers wrote:
> I have a question concerning the new Windows toolchain for R>= 2.14.2.
> When trying out my package 'pracma' on the win-builder development
version
> it will stop with the following error message:
>
>    >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<=
1))
>    >  dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err
= 2e-8
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, 
:
>      non-finite function value
>    Calls: dblquad ...
>           <Anonymous>  ->  f ->  do.call ->  mapply -> 
<Anonymous>  ->  integrate
>    Execution halted
>    ** running examples for arch 'x64' ... ERROR
>    Running examples in 'pracma-Ex.R' failed
>
> This probably means that the following expression got negative for some
> values x, y:
>
>    (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
I think you're right, it's a bug, hopefully easy to fix.  Here's a 
simpler version:

x <- 0*(-1)
sqrt(x)

x is a "negative zero", and the sqrt() function incorrectly produces a
NaN in the new toolchain.

Duncan Murdoch
>
> It appears to be an often used trick in numerical analysis. One advantage
is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>
> The example runs fine on Debian Linux and Mac OS X 32-/64-bit
architectures.
> In my understanding the approach is correct and, as said above, often used
in
> numerical applications.
>
> Can someone explain to me why this fails for the Windows 64-bit compiler
and
> what I should use instead. Thanks.
>
> Hans Werner Borchers
> ABB Corporate Research
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

On Jan 27, 2012, at 13:23 , Hans W Borchers wrote:
>  (1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)
> 
> It appears to be an often used trick in numerical analysis. One advantage
is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
However, "pmax(0, 1 - (x^2 + y^2))" is (unless 0-length x,y is an
issue).

But of course, Duncan is right: It is a bug if you can't take the square
root of negative zero.

-- 
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

This did turn out to be a bug in the new toolchain, and Brian Ripley has 
devised a patch and put together a new one.  I've uploaded a new 
Rtools215.exe, which should be available for download tomorrow, and 
builds of R-patched and R-devel will soon use it.  Everything takes a 
while to propagate to the volunteers and systems that build binaries and 
the mirrors, but we should all be up to date by the end of the week or so.

Thanks for the report!

Duncan Murdoch


On 27/01/2012 7:23 AM, Hans W Borchers wrote:
> I have a question concerning the new Windows toolchain for R>= 2.14.2.
> When trying out my package 'pracma' on the win-builder development
version
> it will stop with the following error message:
>
>    >  f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<=
1))
>    >  dblquad(f3, -1, 1, -1, 1)     #   2.094395124 , i.e. 2/3*pi , err
= 2e-8
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>    Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, 
:
>      non-finite function value
>    Calls: dblquad ...
>           <Anonymous>  ->  f ->  do.call ->  mapply -> 
<Anonymous>  ->  integrate
>    Execution halted
>    ** running examples for arch 'x64' ... ERROR
>    Running examples in 'pracma-Ex.R' failed
>
> This probably means that the following expression got negative for some
> values x, y:
>
>    (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
>
> It appears to be an often used trick in numerical analysis. One advantage
is
> that a function using it is immediately vectorized while an expression such
> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>
> The example runs fine on Debian Linux and Mac OS X 32-/64-bit
architectures.
> In my understanding the approach is correct and, as said above, often used
in
> numerical applications.
>
> Can someone explain to me why this fails for the Windows 64-bit compiler
and
> what I should use instead. Thanks.
>
> Hans Werner Borchers
> ABB Corporate Research
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel