Re: [Rd] Numerical accuracy of matrix multiplication

2016-09-20 Thread Martin Maechler
> peter dalgaard 
> on Fri, 16 Sep 2016 13:33:11 +0200 writes:

> On 16 Sep 2016, at 12:41 , Alexis Sarda  wrote:

>> Hello,
>> 
>> while testing the crossprod() function under Linux, I noticed the 
following:
>> 
>> set.seed(883)
>> x <- rnorm(100)
>> x %*% x - sum(x^2) # equal to 1.421085e-14
>> 
>> Is this difference normal? It seems to be rather large for double 
precision.
>> 

> It's less than .Machine$double.eps, relative (!) to x  %*% x ~= 100.

indeed!

Still, it gives exactly 0 on my platform(s), where I'm using R's
own version of BLAS / Lapack.

Are you perhaps using an "optimized" BLAS / LAPACK , i.e, one
that is fast but slightly less so accurate ?

Martin Maechler,
ETH Zurich


> -pd

>> Regards,
>> Alexis.
>> 
>> [[alternative HTML version deleted]]
>> 
>> __
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel

> -- 
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd@cbs.dk  Priv: pda...@gmail.com

> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


Re: [Rd] Numerical accuracy of matrix multiplication

2016-09-20 Thread Martin Maechler
> Alexis Sarda 
> on Tue, 20 Sep 2016 17:33:49 +0200 writes:

> I just realized that I was actually using a different random number
> generator, could that be a valid reason for the discrepancy?

> The code should be:

> RNGkind("L'Ecuyer")
> set.seed(883)
> x <- rnorm(100)
> x %*% x - sum(x^2) # equal to 1.421085e-14


Yes, now I get the same result so my story on  "BLAS / LAPACK"
is not relevant here.

But do note the main point from Peter Dalgaard that this is well
within Machine epsilon precision.

More precisely, here it is really one bit difference in the
least significant bit :

> print(rbind( x%*%x, crossprod(x), sum(x^2)), digits= 19)
 [,1]
[1,] 103.5096830356289814
[2,] 103.5096830356289814
[3,] 103.5096830356289672
> cbind(sprintf("%a", c(x%*%x, crossprod(x), sum(x^2
 [,1]  
[1,] "0x1.9e09ea598568fp+6"
[2,] "0x1.9e09ea598568fp+6"
[3,] "0x1.9e09ea598568ep+6"
> 


> Regards,
> Alexis Sarda.



> On Tue, Sep 20, 2016 at 5:27 PM, Martin Maechler 
> wrote:

>> > peter dalgaard 
>> > on Fri, 16 Sep 2016 13:33:11 +0200 writes:
>> 
>> > On 16 Sep 2016, at 12:41 , Alexis Sarda 
>> wrote:
>> 
>> >> Hello,
>> >>
>> >> while testing the crossprod() function under Linux, I noticed the
>> following:
>> >>
>> >> set.seed(883)
>> >> x <- rnorm(100)
>> >> x %*% x - sum(x^2) # equal to 1.421085e-14
>> >>
>> >> Is this difference normal? It seems to be rather large for double
>> precision.
>> >>
>> 
>> > It's less than .Machine$double.eps, relative (!) to x  %*% x ~= 100.
>> 
>> indeed!
>> 
>> Still, it gives exactly 0 on my platform(s), where I'm using R's
>> own version of BLAS / Lapack.
>> 
>> Are you perhaps using an "optimized" BLAS / LAPACK , i.e, one
>> that is fast but slightly less so accurate ?
>> 
>> Martin Maechler,
>> ETH Zurich
>> 
>> 
>> > -pd
>> 
>> >> Regards,
>> >> Alexis.
>> >>
>> >> [[alternative HTML version deleted]]
>> >>
>> >> __
>> >> R-devel@r-project.org mailing list
>> >> https://stat.ethz.ch/mailman/listinfo/r-devel
>> 
>> > --
>> > Peter Dalgaard, Professor,
>> > Center for Statistics, Copenhagen Business School
>> > Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>> > Phone: (+45)38153501
>> > Office: A 4.23
>> > Email: pd@cbs.dk  Priv: pda...@gmail.com
>> 
>> > __
>> > R-devel@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-devel
>> 

> [[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


Re: [Rd] Numerical accuracy of matrix multiplication

2016-09-20 Thread Alexis Sarda
I just realized that I was actually using a different random number
generator, could that be a valid reason for the discrepancy?

The code should be:

RNGkind("L'Ecuyer")
set.seed(883)
x <- rnorm(100)
x %*% x - sum(x^2) # equal to 1.421085e-14


Regards,
Alexis Sarda.



On Tue, Sep 20, 2016 at 5:27 PM, Martin Maechler  wrote:

> > peter dalgaard 
> > on Fri, 16 Sep 2016 13:33:11 +0200 writes:
>
> > On 16 Sep 2016, at 12:41 , Alexis Sarda 
> wrote:
>
> >> Hello,
> >>
> >> while testing the crossprod() function under Linux, I noticed the
> following:
> >>
> >> set.seed(883)
> >> x <- rnorm(100)
> >> x %*% x - sum(x^2) # equal to 1.421085e-14
> >>
> >> Is this difference normal? It seems to be rather large for double
> precision.
> >>
>
> > It's less than .Machine$double.eps, relative (!) to x  %*% x ~= 100.
>
> indeed!
>
> Still, it gives exactly 0 on my platform(s), where I'm using R's
> own version of BLAS / Lapack.
>
> Are you perhaps using an "optimized" BLAS / LAPACK , i.e, one
> that is fast but slightly less so accurate ?
>
> Martin Maechler,
> ETH Zurich
>
>
> > -pd
>
> >> Regards,
> >> Alexis.
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> __
> >> R-devel@r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> > --
> > Peter Dalgaard, Professor,
> > Center for Statistics, Copenhagen Business School
> > Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> > Phone: (+45)38153501
> > Office: A 4.23
> > Email: pd@cbs.dk  Priv: pda...@gmail.com
>
> > __
> > R-devel@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
>

[[alternative HTML version deleted]]

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


Re: [Rd] Handlers in setGraphicsEventHandlers() can recursively call getGraphicsEvent(). Intended behavior?

2016-09-20 Thread Paul Murrell

Hi

Is the correct patch to remove the setting of the gettingEvent flag or 
would it be better to flip the TRUE/FALSE setting (set to TRUE before 
handling then reset to FALSE after handling) ?


Also, for this patch and for the other two you sent, one difficulty will 
be with testing the patches.  I have no testing code for this, so would 
need at least a test or two from you (ideally someone would also have 
some regression tests, beyond ?getGraphicsEvent, to ensure continuity of 
previous behaviour).


Paul

On 18/09/2016 3:29 a.m., Richard Bodewits wrote:

Hey all.

As in general it's a bad idea to allow an event handler to generate an
event, and as comments in the code seem to suggest this isn't the
intention either, I was wondering about recursion in getGraphicsEvent().
In main/gevents.c, both doMouseEvent() and doKeybd() have the following
line;

dd->gettingEvent = FALSE; /* avoid recursive calls */

And they reset it to TRUE before returning. The effective result of this
is that the event handlers on the R side are allowed to call
getGraphicsEvent(), and recurse until they eventually would run out of
stack space. Though R does catch and handle that cleanly with the error;

Error: evaluation nested too deeply: infinite recursion /
options(expressions=)?

A quick scan of the SVN logs suggests this code has been untouched since
its introduction in 2004, so I'm left to wonder if this is intended
behavior. It stands out as a bit of a sore thumb due to the generic
check for recursion in do_getGraphicsEvent() in the same file, which
will error out with error(_("recursive use of 'getGraphicsEvent' not
supported")) if dd->gettingEvent is already set to TRUE. Which would
suggest recursively calling it is very much not intended to be possible.

To me, setting gettingEvent to FALSE seems like an easy mistake to make
if you temporarily interpret gettingEvent to mean that event(s) are
allowed to still come in. But the actual interpretation in
do_getGraphicsEvents() is the opposite, as it's interpreted as an
indicator of whether or not an event is currently being processed.

I've removed the gettingEvent altering lines from both doMouseEvent()
and doKeybd() to no ill effect, and doing so disabled the ability to
call getGraphicsEvent() from inside one of the assigned handlers as
expected. But after 12 (!) years, it's conceivable that people have come
to depend on this behavior in existing scripts. Is this something that
should be left alone to minimize disruption? Or should this be fixed (if
it is indeed unintended) for the sake of protecting people from infinite
recursions?

I've included a small patch as attachment that removes the offending lines.


- Richard Bodewits



__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel



--
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


[Rd] Undocumented 'use.names' argument to c()

2016-09-20 Thread Karl Millar via R-devel
'c' has an undocumented 'use.names' argument.  I'm not sure if this is
a documentation or implementation bug.

> c(a = 1)
a
1
> c(a = 1, use.names = F)
[1] 1

Karl

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


Re: [Rd] Undocumented 'use.names' argument to c()

2016-09-20 Thread David Winsemius

> On Sep 20, 2016, at 7:18 PM, Karl Millar via R-devel  
> wrote:
> 
> 'c' has an undocumented 'use.names' argument.  I'm not sure if this is
> a documentation or implementation bug.

It came up on stackoverflow a couple of years ago:

http://stackoverflow.com/questions/24815572/why-does-function-c-accept-an-undocumented-argument/24815653#24815653

At the time it appeared to me to be a documentation lag.


> 
>> c(a = 1)
> a
> 1
>> c(a = 1, use.names = F)
> [1] 1
> 
> Karl
> 
> __
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

David Winsemius
Alameda, CA, USA

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel