G'day Martin,
On Mon, 3 Mar 2008 10:16:45 +0100
Martin Maechler <[EMAIL PROTECTED]> wrote:
> >>>>> "BAT" == Berwin A Turlach <[EMAIL PROTECTED]>
> >>>>> on Fri, 29 Feb 2008 18:19:40 +0800 writes:
>
> BAT> while looking for some inspiration of how to organise some
> BAT> code, I studied the code of random.c and noticed that for
> BAT> distributions with 2 or 3 parameters the user is not warned
> BAT> if NAs are created while such a warning is issued for
> BAT> distributions with 1 parameter. [...] The attached patch
> BAT> rectifies this. [...]
>
> I cannot imagine a design reason for that. If there was, it
> should have been mentioned as a comment in the C code.
>
> I'll commit your patch (if it passes the checks).
Sorry, I was a bit in a hurry when writing the e-mail, so I forgot to
mention that the source code modified by this patch compiles and passes
"make check FORCE=FORCE" on my machine.
And in my hurry, I also posted from my NUS account, without realising
it, which forced you to intervene as moderator and to approve the
posting. My apologies for the extra work. But this gave me the idea
to also subscribe to r-devel with my NUS account and configure the
subscriptions so that I only receive e-mail at my UWA account. Thus,
hopefully, you will not have to intervene again. (Which this e-mail
should test.)
> BAT> BTW, there are other places in the code were NAs can be
> BAT> created but no warning is issued. E.g:
>
> >> rexp(2, rate=numeric())
> BAT> [1] NA NA
> >> rnorm(2, mean=numeric())
> BAT> [1] NA NA
>
> BAT> I wonder whether a warning should be issued in this case
> BAT> too.
>
> Yes, "should in principle".
>
> If you feel like finding another elegant patch...
Well, elegance is in the eye of the beholder. :-) I attach two
patches. One that adds warning messages at the other places where NAs
can be generated.
The second one in additiona rearranges the code a bit such that in the
case when all the "vectors" that contain the parameter values of the
distribution, from which one wants to simulate, are of length one some
unnecessary calculations is taken out of the for loop. I am not sure
how much time is actually saved in this situation, but I belong to the
school that things such kind of optimisation should be done. :) If you
think it bloats the code too much (or duplicates things too much
leading to hard to maintain code), then feel free to ignore this second
patch.
> Thank you Berwin, for your contribution!
My pleasure.
Cheers,
Berwin
Index: src/main/random.c
===================================================================
--- src/main/random.c (revision 44677)
+++ src/main/random.c (working copy)
@@ -81,6 +81,7 @@
if (na < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(CADR(args), REALSXP));
@@ -123,7 +124,7 @@
#define RAND2(num,name) \
case num: \
- random2(name, REAL(a), na, REAL(b), nb, REAL(x), n); \
+ naflag = random2(name, REAL(a), na, REAL(b), nb, REAL(x), n); \
break
/* "do_random2" - random sampling from 2 parameter families. */
@@ -155,6 +156,7 @@
if (na < 1 || nb < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(CADR(args), REALSXP));
@@ -207,7 +209,7 @@
#define RAND3(num,name) \
case num: \
- random3(name, REAL(a), na, REAL(b), nb, REAL(c), nc, REAL(x),
n); \
+ naflag = random3(name, REAL(a), na, REAL(b), nb, REAL(c), nc,
REAL(x), n); \
break
@@ -244,6 +246,7 @@
if (na < 1 || nb < 1 || nc < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(a, REALSXP));
Index: src/main/random.c
===================================================================
--- src/main/random.c (revision 44677)
+++ src/main/random.c (working copy)
@@ -41,10 +41,19 @@
double ai;
int i;
errno = 0;
- for (i = 0; i < n; i++) {
+ if( na == 1 ){
+ ai = a[0];
+ for (i = 0; i < n; i++) {
+ x[i] = f(ai);
+ if (!R_FINITE(x[i])) naflag = 1;
+ }
+ }
+ else{
+ for (i = 0; i < n; i++) {
ai = a[i % na];
x[i] = f(ai);
if (!R_FINITE(x[i])) naflag = 1;
+ }
}
return(naflag);
}
@@ -81,6 +90,7 @@
if (na < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(CADR(args), REALSXP));
@@ -112,18 +122,28 @@
double ai, bi; int i;
Rboolean naflag = FALSE;
errno = 0;
- for (i = 0; i < n; i++) {
+ if( na == 1 && na == 1 ){
+ ai = a[0];
+ bi = b[0];
+ for (i = 0; i < n; i++) {
+ x[i] = f(ai, bi);
+ if (!R_FINITE(x[i])) naflag = 1;
+ }
+ }
+ else{
+ for (i = 0; i < n; i++) {
ai = a[i % na];
bi = b[i % nb];
x[i] = f(ai, bi);
if (!R_FINITE(x[i])) naflag = 1;
+ }
}
return(naflag);
}
#define RAND2(num,name) \
case num: \
- random2(name, REAL(a), na, REAL(b), nb, REAL(x), n); \
+ naflag = random2(name, REAL(a), na, REAL(b), nb, REAL(x), n); \
break
/* "do_random2" - random sampling from 2 parameter families. */
@@ -155,6 +175,7 @@
if (na < 1 || nb < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(CADR(args), REALSXP));
@@ -195,19 +216,30 @@
int i;
Rboolean naflag = FALSE;
errno = 0;
- for (i = 0; i < n; i++) {
+ if( na == 1 && nb == 1 && nc == 1 ){
+ ai = a[0];
+ bi = b[0];
+ ci = c[0];
+ for (i = 0; i < n; i++) {
+ x[i] = f(ai, bi, ci);
+ if (!R_FINITE(x[i])) naflag = TRUE;
+ }
+ }
+ else{
+ for (i = 0; i < n; i++) {
ai = a[i % na];
bi = b[i % nb];
ci = c[i % nc];
x[i] = f(ai, bi, ci);
if (!R_FINITE(x[i])) naflag = TRUE;
+ }
}
return(naflag);
}
#define RAND3(num,name) \
case num: \
- random3(name, REAL(a), na, REAL(b), nb, REAL(c), nc, REAL(x),
n); \
+ naflag = random3(name, REAL(a), na, REAL(b), nb, REAL(c), nc,
REAL(x), n); \
break
@@ -244,6 +276,7 @@
if (na < 1 || nb < 1 || nc < 1) {
for (i = 0; i < n; i++)
REAL(x)[i] = NA_REAL;
+ warning(_("NAs produced"));
}
else {
PROTECT(a = coerceVector(a, REALSXP));
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
