[Pomp-commits] r860 - in pkg/pomp: . R inst inst/include src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 7 14:13:08 CEST 2013
Author: kingaa
Date: 2013-06-07 14:13:08 +0200 (Fri, 07 Jun 2013)
New Revision: 860
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/include/pomp.h
pkg/pomp/src/R_init_pomp.c
pkg/pomp/src/eulermultinom.c
pkg/pomp/src/lookup_table.c
pkg/pomp/src/pomp.h
Log:
- reulermultinom, deulermultinom, and dot_product are no longer linkables exported in src/R_init_pomp.c but rather are defined as static inlines in the pomp.h header file.
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/DESCRIPTION 2013-06-07 12:13:08 UTC (rev 860)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.45-3
-Date: 2013-06-05
+Version: 0.45-4
+Date: 2013-06-06
Maintainer: Aaron A. King <kingaa at umich.edu>
Authors at R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/R/builder.R 2013-06-07 12:13:08 UTC (rev 860)
@@ -107,7 +107,6 @@
decl <- list(
periodic_bspline_basis_eval="\tvoid (*periodic_bspline_basis_eval)(double,double,int,int,double*);\nperiodic_bspline_basis_eval = (void (*)(double,double,int,int,double*)) R_GetCCallable(\"pomp\",\"periodic_bspline_basis_eval\");\n",
- reulermultinom="\tvoid (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n",
get_pomp_userdata="\tconst SEXP (*get_pomp_userdata)(const char *);\npomp_get_userdata = (const SEXP (*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata\");\n",
get_pomp_userdata_int="\tconst int * (*get_pomp_userdata_int)(const char *);\npomp_get_userdata_int = (const int *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_int\");\n",
get_pomp_userdata_double="\tconst double * (*get_pomp_userdata_double)(const char *);\npomp_get_userdata_double = (const double *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_double\");\n"
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/inst/NEWS 2013-06-07 12:13:08 UTC (rev 860)
@@ -1,4 +1,7 @@
NEWS
+0.45-4
+ o changes in the way 'reulermultinom', 'deulermultinom', 'dot_product' are exported to other packages. Rather than being exported as linkables, these are now defined as static inline functions in the 'pomp.h' header.
+
0.45-3
o fix bug with 'continue' and method 'mif2'
Modified: pkg/pomp/inst/include/pomp.h
===================================================================
--- pkg/pomp/inst/include/pomp.h 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/inst/include/pomp.h 2013-06-07 12:13:08 UTC (rev 860)
@@ -5,15 +5,15 @@
#include <R.h>
#include <Rmath.h>
-#include <Rdefines.h>
+#include <Rinternals.h>
#include <R_ext/Rdynload.h>
-// facility for extracting R objects from the 'userdata' slot
+// facilities for extracting R objects from the 'userdata' slot
const SEXP get_pomp_userdata (const char *name);
const int *get_pomp_userdata_int (const char *name);
const double *get_pomp_userdata_double (const char *name);
-// facility for computing evaluating a basis of periodic bsplines
+// facility for evaluating a set of periodic bspline basis functions
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
// Prototype for parameter transformation function.
@@ -191,7 +191,13 @@
// facility for computing the inner product of
// a vector of parameters ('coef') against a vector of basis-function values ('basis')
-double dot_product (int dim, const double *basis, const double *coef);
+static R_INLINE double dot_product (int dim, const double *basis, const double *coef) {
+ int j;
+ double trans = 0.0;
+ for (j = 0; j < dim; j++)
+ trans += coef[j]*basis[j];
+ return(trans);
+}
static R_INLINE double logit (double p) {
return log(p/(1.0-p));
@@ -201,17 +207,88 @@
return 1.0/(1.0+exp(-x));
}
-// prototypes for C-level access to Euler-multinomial distribution functions
+// C-level definitions of Euler-multinomial distribution functions
// simulate Euler-multinomial transitions
// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally
// this must be done by the calling program
// But note that when reulermultinom is called inside a pomp 'rprocess', there is no need to call
// {Get,Put}RNGState() as this is handled by pomp
-void reulermultinom (int m, double size, double *rate, double dt, double *trans);
+static void reulermultinom (int m, double size, double *rate, double dt, double *trans) {
+ double p = 0.0;
+ int j, k;
+ if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
+ for (k = 0; k < m; k++) trans[k] = R_NaN;
+ return;
+ }
+ for (k = 0; k < m; k++) {
+ if (rate[k] < 0.0) {
+ for (j = 0; j < m; j++) trans[j] = R_NaN;
+ return;
+ }
+ p += rate[k]; // total event rate
+ }
+ if (p > 0.0) {
+ size = rbinom(size,1-exp(-p*dt)); // total number of events
+ if (!(R_FINITE(size)))
+ warning("reulermultinom: result of binomial draw is not finite");
+ m -= 1;
+ for (k = 0; k < m; k++) {
+ if (rate[k] > p) p = rate[k];
+ trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
+ if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
+ warning("reulermultinom: result of binomial draw is not finite");
+ size -= trans[k];
+ p -= rate[k];
+ }
+ trans[m] = size;
+ } else {
+ for (k = 0; k < m; k++) trans[k] = 0.0;
+ }
+}
-// compute probabilities of eulermultinomial transitions
-double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log);
+// compute probabilities of Euler-multinomial transitions
+static double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log) {
+ double p = 0.0;
+ double n = 0.0;
+ double ff = 0.0;
+ int k;
+ if ((dt < 0.0) || (size < 0.0) || (floor(size+0.5) != size)) {
+ warning("NaNs produced");
+ return R_NaN;
+ }
+ for (k = 0; k < m; k++) {
+ if (rate[k] < 0.0) {
+ warning("NaNs produced");
+ return R_NaN;
+ }
+ if (trans[k] < 0.0) {
+ ff = (give_log) ? R_NegInf: 0.0;
+ return ff;
+ }
+ p += rate[k]; // total event rate
+ n += trans[k]; // total number of events
+ }
+ if (n > size) {
+ ff = (give_log) ? R_NegInf: 0.0;
+ return ff;
+ }
+ ff = dbinom(n,size,1-exp(-p*dt),1); // total number of events
+ m -= 1;
+ for (k = 0; k < m; k++) {
+ if ((n > 0) && (p > 0)) {
+ if (rate[k] > p) p = rate[k];
+ ff += dbinom(trans[k],n,rate[k]/p,1);
+ } else if (trans[k] > 0.0) {
+ ff = R_NegInf;
+ return ff;
+ }
+ n -= trans[k];
+ p -= rate[k];
+ }
+ ff = (give_log) ? ff : exp(ff);
+ return ff;
+}
static R_INLINE double rbetabinom (double size, double prob, double theta) {
return rbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
Modified: pkg/pomp/src/R_init_pomp.c
===================================================================
--- pkg/pomp/src/R_init_pomp.c 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/src/R_init_pomp.c 2013-06-07 12:13:08 UTC (rev 860)
@@ -4,9 +4,6 @@
void R_init_pomp (DllInfo *info) {
R_RegisterCCallable("pomp","periodic_bspline_basis_eval",(DL_FUNC) &periodic_bspline_basis_eval);
- R_RegisterCCallable("pomp","dot_product",(DL_FUNC) &dot_product);
- R_RegisterCCallable("pomp","reulermultinom",(DL_FUNC) &reulermultinom);
- R_RegisterCCallable("pomp","deulermultinom",(DL_FUNC) &deulermultinom);
R_RegisterCCallable("pomp","get_pomp_userdata",(DL_FUNC) &get_pomp_userdata);
R_RegisterCCallable("pomp","get_pomp_userdata_int",(DL_FUNC) &get_pomp_userdata_int);
R_RegisterCCallable("pomp","get_pomp_userdata_double",(DL_FUNC) &get_pomp_userdata_double);
Modified: pkg/pomp/src/eulermultinom.c
===================================================================
--- pkg/pomp/src/eulermultinom.c 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/src/eulermultinom.c 2013-06-07 12:13:08 UTC (rev 860)
@@ -4,82 +4,6 @@
#include "pomp_internal.h"
-void reulermultinom (int m, double size, double *rate, double dt, double *trans) {
- double p = 0.0;
- int j, k;
- if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
- for (k = 0; k < m; k++) trans[k] = R_NaN;
- return;
- }
- for (k = 0; k < m; k++) {
- if (rate[k] < 0.0) {
- for (j = 0; j < m; j++) trans[j] = R_NaN;
- return;
- }
- p += rate[k]; // total event rate
- }
- if (p > 0.0) {
- size = rbinom(size,1-exp(-p*dt)); // total number of events
- if (!(R_FINITE(size)))
- warning("reulermultinom: result of binomial draw is not finite");
- m -= 1;
- for (k = 0; k < m; k++) {
- if (rate[k] > p) p = rate[k];
- trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
- if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
- warning("reulermultinom: result of binomial draw is not finite");
- size -= trans[k];
- p -= rate[k];
- }
- trans[m] = size;
- } else {
- for (k = 0; k < m; k++) trans[k] = 0.0;
- }
-}
-
-// probability density of Euler-multinomial transitions
-double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log) {
- double p = 0.0;
- double n = 0.0;
- double ff = 0.0;
- int k;
- if ((dt < 0.0) || (size < 0.0) || (floor(size+0.5) != size)) {
- warning("NaNs produced");
- return R_NaN;
- }
- for (k = 0; k < m; k++) {
- if (rate[k] < 0.0) {
- warning("NaNs produced");
- return R_NaN;
- }
- if (trans[k] < 0.0) {
- ff = (give_log) ? R_NegInf: 0.0;
- return ff;
- }
- p += rate[k]; // total event rate
- n += trans[k]; // total number of events
- }
- if (n > size) {
- ff = (give_log) ? R_NegInf: 0.0;
- return ff;
- }
- ff = dbinom(n,size,1-exp(-p*dt),1); // total number of events
- m -= 1;
- for (k = 0; k < m; k++) {
- if ((n > 0) && (p > 0)) {
- if (rate[k] > p) p = rate[k];
- ff += dbinom(trans[k],n,rate[k]/p,1);
- } else if (trans[k] > 0.0) {
- ff = R_NegInf;
- return ff;
- }
- n -= trans[k];
- p -= rate[k];
- }
- ff = (give_log) ? ff : exp(ff);
- return ff;
-}
-
void reulermultinom_multi (int *n, int *ntrans, double *size, double *rate, double *dt, double *trans) {
int k;
int m = *ntrans;
Modified: pkg/pomp/src/lookup_table.c
===================================================================
--- pkg/pomp/src/lookup_table.c 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/src/lookup_table.c 2013-06-07 12:13:08 UTC (rev 860)
@@ -70,13 +70,3 @@
}
}
-// facility for computing the inner product of
-// a vector of parameters ('coef') against a vector of basis-function values ('basis')
-double dot_product (int dim, const double *basis, const double *coef) {
- int j;
- double trans = 0.0;
- for (j = 0; j < dim; j++)
- trans += coef[j]*basis[j];
- return(trans);
-}
-
Modified: pkg/pomp/src/pomp.h
===================================================================
--- pkg/pomp/src/pomp.h 2013-06-06 00:03:28 UTC (rev 859)
+++ pkg/pomp/src/pomp.h 2013-06-07 12:13:08 UTC (rev 860)
@@ -5,15 +5,15 @@
#include <R.h>
#include <Rmath.h>
-#include <Rdefines.h>
+#include <Rinternals.h>
#include <R_ext/Rdynload.h>
-// facility for extracting R objects from the 'userdata' slot
+// facilities for extracting R objects from the 'userdata' slot
const SEXP get_pomp_userdata (const char *name);
const int *get_pomp_userdata_int (const char *name);
const double *get_pomp_userdata_double (const char *name);
-// facility for computing evaluating a basis of periodic bsplines
+// facility for evaluating a set of periodic bspline basis functions
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
// Prototype for parameter transformation function.
@@ -191,7 +191,13 @@
// facility for computing the inner product of
// a vector of parameters ('coef') against a vector of basis-function values ('basis')
-double dot_product (int dim, const double *basis, const double *coef);
+static R_INLINE double dot_product (int dim, const double *basis, const double *coef) {
+ int j;
+ double trans = 0.0;
+ for (j = 0; j < dim; j++)
+ trans += coef[j]*basis[j];
+ return(trans);
+}
static R_INLINE double logit (double p) {
return log(p/(1.0-p));
@@ -201,17 +207,88 @@
return 1.0/(1.0+exp(-x));
}
-// prototypes for C-level access to Euler-multinomial distribution functions
+// C-level definitions of Euler-multinomial distribution functions
// simulate Euler-multinomial transitions
// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally
// this must be done by the calling program
// But note that when reulermultinom is called inside a pomp 'rprocess', there is no need to call
// {Get,Put}RNGState() as this is handled by pomp
-void reulermultinom (int m, double size, double *rate, double dt, double *trans);
+static void reulermultinom (int m, double size, double *rate, double dt, double *trans) {
+ double p = 0.0;
+ int j, k;
+ if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
+ for (k = 0; k < m; k++) trans[k] = R_NaN;
+ return;
+ }
+ for (k = 0; k < m; k++) {
+ if (rate[k] < 0.0) {
+ for (j = 0; j < m; j++) trans[j] = R_NaN;
+ return;
+ }
+ p += rate[k]; // total event rate
+ }
+ if (p > 0.0) {
+ size = rbinom(size,1-exp(-p*dt)); // total number of events
+ if (!(R_FINITE(size)))
+ warning("reulermultinom: result of binomial draw is not finite");
+ m -= 1;
+ for (k = 0; k < m; k++) {
+ if (rate[k] > p) p = rate[k];
+ trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
+ if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
+ warning("reulermultinom: result of binomial draw is not finite");
+ size -= trans[k];
+ p -= rate[k];
+ }
+ trans[m] = size;
+ } else {
+ for (k = 0; k < m; k++) trans[k] = 0.0;
+ }
+}
-// compute probabilities of eulermultinomial transitions
-double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log);
+// compute probabilities of Euler-multinomial transitions
+static double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log) {
+ double p = 0.0;
+ double n = 0.0;
+ double ff = 0.0;
+ int k;
+ if ((dt < 0.0) || (size < 0.0) || (floor(size+0.5) != size)) {
+ warning("NaNs produced");
+ return R_NaN;
+ }
+ for (k = 0; k < m; k++) {
+ if (rate[k] < 0.0) {
+ warning("NaNs produced");
+ return R_NaN;
+ }
+ if (trans[k] < 0.0) {
+ ff = (give_log) ? R_NegInf: 0.0;
+ return ff;
+ }
+ p += rate[k]; // total event rate
+ n += trans[k]; // total number of events
+ }
+ if (n > size) {
+ ff = (give_log) ? R_NegInf: 0.0;
+ return ff;
+ }
+ ff = dbinom(n,size,1-exp(-p*dt),1); // total number of events
+ m -= 1;
+ for (k = 0; k < m; k++) {
+ if ((n > 0) && (p > 0)) {
+ if (rate[k] > p) p = rate[k];
+ ff += dbinom(trans[k],n,rate[k]/p,1);
+ } else if (trans[k] > 0.0) {
+ ff = R_NegInf;
+ return ff;
+ }
+ n -= trans[k];
+ p -= rate[k];
+ }
+ ff = (give_log) ? ff : exp(ff);
+ return ff;
+}
static R_INLINE double rbetabinom (double size, double prob, double theta) {
return rbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
More information about the pomp-commits
mailing list