Commit c4aad1f4 authored by Radford Neal's avatar Radford Neal

work on gradients of random functions, bug fixes relating to gradients

parent af980c5f
......@@ -497,7 +497,7 @@ typedef struct {
setprim_ptr->primsxp.whole \
= (R_FunTab[setprim_value].eval/1000000)&1; \
setprim_ptr->primsxp.gradn \
= (R_FunTab[setprim_value].eval/10000000)&3; \
= (R_FunTab[setprim_value].eval/10000000)%10; \
} while (0)
#define PRIMOFFSET(x) \
......
......@@ -376,9 +376,11 @@ struct primsxp_struct { /* table offset of this and other info is in gp */
unsigned int fast_sub:1; /* subassign/subset fn that can use fast method*/
unsigned int dsptch1:1; /* might dispatch on 1st argument (only
for when fast_cfun != NULL) */
unsigned int gradn:2; /* # of internal args that might have gradient */
unsigned int whole:1; /* Do special processing for .Internal
when VARIANT_WHOLE_BODY (BUILTIN only) */
unsigned int padbit:1;
unsigned int gradn:3; /* # of internal args that might have gradient,
(with fudge, see 'eval' field in names.c) */
#if USE_COMPRESSED_POINTERS && SIZEOF_CHAR_P == 4
int32_t padding1, padding2;
#endif
......
......@@ -2807,9 +2807,9 @@ SEXP attribute_hidden evalList_gradient (SEXP el, SEXP rho, int variant, int n)
int i, m;
i = 0;
if (n > 4) {
m = 0; /* so i won't be less than m for first argument */
n -= 4; /* m will be set to this new n after first argument */
if (n > 4) { /* don't ask for gradient for first argument */
m = 0; /* so i won't be less than m for first argument */
n -= 3; /* m will be set to this new n after first argument */
}
else
m = n;
......@@ -2852,7 +2852,7 @@ SEXP attribute_hidden evalList_gradient (SEXP el, SEXP rho, int variant, int n)
} else {
if (CDR(el) == R_NilValue)
varpend = variant; /* don't defer pointlessly for last one */
varpend = variant; /* don't defer pointlessly for last one */
INC_NAMEDCNT(CAR(tail)); /* OK when tail is R_NilValue */
ev_el = EVALV (CAR(el), rho,
i<m ? varpend | VARIANT_GRADIENT : varpend);
......
......@@ -1629,7 +1629,8 @@ SEXP attribute_hidden mkValuePROMISE(SEXP expr, SEXP value)
}
/* mkPRIMSXP - return a primitve function, "builtin" or "special".
/* mkPRIMSXP - Return a primitve function, "builtin" or "special".
May be actually "primitive", or be "internal".
Primitive objects are recorded to avoid creation of extra ones
during unserializaton or reconstruction after a package has
......
......@@ -598,14 +598,15 @@ static void SetupBuiltins(void)
int this_code = R_FunTab[i].code;
SEXP prim;
/* prim needs protect since install can (and does here) allocate.
Except... mkPRIMSXP now caches them all, so maybe not. */
Except... they're now uncollected... */
PROTECT(prim = mkPRIMSXP(i, R_FunTab[i].eval % 10));
if (0) { /* enable to display info */
Rprintf (
"SETUP: %s, builtin %d, internal %d, visible %d, pending_ok %d\n",
PRIMNAME(prim), TYPEOF(prim) == BUILTINSXP,
PRIMINTERNAL(prim), PRIMPRINT(prim), PRIMPENDINGOK(prim));
"SETUP: %d %s, builtin %d, int %d, vis %d, pend_ok %d, gradn %d\n",
i, PRIMNAME(prim), TYPEOF(prim) == BUILTINSXP,
PRIMINTERNAL(prim), PRIMPRINT(prim), PRIMPENDINGOK(prim),
PRIMGRADN(prim));
}
if (TYPEOF(prim) == SPECIALSXP && !PRIMINTERNAL(prim)
......
......@@ -41,10 +41,25 @@ R_NORETURN static void invalid(SEXP call)
error(_("invalid arguments"));
}
/* "do_random1" - random sampling from 1 parameter families. */
/* Derivatives for random1 generators. */
static double (*rand1_funs[6])(double) = {
rchisq, rexp, rgeom, rpois, rt, rsignrank
static double Drexp (double r, double scale)
{
return r / scale;
}
/* Table of functions for generation and derivatives of random1 distibutions. */
static struct {
double (*fncall)(double); double (*Dcall)(double,double);
} rand1_table[6] =
{
{ rchisq, 0 },
{ rexp, Drexp },
{ rgeom, 0 },
{ rpois, 0 },
{ rt, 0 },
{ rsignrank,0 }
};
static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
......@@ -70,15 +85,34 @@ static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
na = LENGTH(a);
int opcode = PRIMVAL(op);
if (opcode < 0 || opcode >= sizeof rand1_funs)
if (opcode < 0 || opcode > 5)
error(_("internal error in do_random1"));
double (*rf)(double) = rand1_funs[opcode];
double (*fncall)(double) = rand1_table[opcode].fncall;
if (n == 1 && na >= 1) { /* quickly generate single value */
GetRNGstate();
double r = rf (*REAL(a));
if (ISNAN(r)) warning(_("NAs produced"));
double av = *REAL(a);
double r = fncall (av);
PutRNGstate();
if (ISNAN(r)) {
warning(_("NAs produced"));
}
else {
/* Compute gradient if requested. */
SEXP g = ATTRIB(CDR(args));
if (g != R_NilValue) {
double (*Dcall)(double,double) = rand1_table[opcode].Dcall;
if (Dcall != 0) {
R_gradient = copy_scaled_gradients (g, Dcall(r,av));
R_variant_result = VARIANT_GRADIENT_FLAG;
}
}
}
UNPROTECT(1); /* a */
return ScalarReal(r);
}
......@@ -107,10 +141,10 @@ static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
double *ap = REAL(a), *xp = REAL(x);
if (na == 1) {
double ar = *ap;
double av = *ap;
int i;
for (i = 0; i < n; i++) {
xp[i] = rf (ar);
xp[i] = fncall (av);
if (ISNAN(xp[i])) naflag = TRUE;
}
}
......@@ -118,7 +152,7 @@ static SEXP do_random1(SEXP call, SEXP op, SEXP args, SEXP rho)
int i, i1;
for (i = 0, i1 = 0; i < n; i++, i1++) {
if (i1 == na) i1 = 0;
xp[i] = rf (ap[i1]);
xp[i] = fncall (ap[i1]);
if (ISNAN(xp[i])) naflag = TRUE;
}
}
......@@ -807,11 +841,11 @@ attribute_hidden FUNTAB R_FunTab_random[] =
{
/* printname c-entry offset eval arity pp-kind precedence rightassoc */
{"rchisq", do_random1, 0, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rexp", do_random1, 1, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rgeom", do_random1, 2, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rpois", do_random1, 3, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rt", do_random1, 4, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rchisq", do_random1, 0, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rexp", do_random1, 1, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rgeom", do_random1, 2, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rpois", do_random1, 3, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rt", do_random1, 4, 51000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rsignrank", do_random1, 5, 1000011, 2, {PP_FUNCALL, PREC_FN, 0}},
{"rbeta", do_random2, 0, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
{"rbinom", do_random2, 1, 1000011, 3, {PP_FUNCALL, PREC_FN, 0}},
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment