[Pomp-commits] r679 - in pkg/pomp: . R inst inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 24 19:09:37 CEST 2012
Author: kingaa
Date: 2012-04-24 19:09:37 +0200 (Tue, 24 Apr 2012)
New Revision: 679
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/inst/ChangeLog
pkg/pomp/inst/NEWS
pkg/pomp/inst/TODO
pkg/pomp/inst/doc/Makefile
pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw
pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
Log:
- put 'pompBuilder' demonstration into 'advanced topics' vignette
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/DESCRIPTION 2012-04-24 17:09:37 UTC (rev 679)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-7
-Date: 2012-04-23
+Version: 0.41-8
+Date: 2012-04-24
Revision: $Rev$
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/R/builder.R 2012-04-24 17:09:37 UTC (rev 679)
@@ -58,12 +58,12 @@
}
pompLink <- function (name) {
- solib <- paste(name,.Platform$dynlib.ext,sep="")
+ solib <- paste0(name,.Platform$dynlib.ext)
dyn.load(solib)
}
pompUnlink <- function (name) {
- solib <- paste(name,.Platform$dynlib.ext,sep="")
+ solib <- paste0(name,.Platform$dynlib.ext)
dyn.unload(solib)
}
@@ -79,7 +79,7 @@
header <- list(
file="/* pomp model file: {%name%} */\n\n#include <pomp.h>\n#include <R_ext/Rdynload.h>\n",
rmeasure="\nvoid {%name%}_rmeasure (double *__y, double *__x, double *__p, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n",
- dmeasure= "\nvoid {%name%}_dmeasure (double *__lik, double *__y, double *__x, double *__p, int __give_log, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n",
+ dmeasure= "\nvoid {%name%}_dmeasure (double *__lik, double *__y, double *__x, double *__p, int give_log, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n",
step.fn="\nvoid {%name%}_stepfn (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, int __covdim, const double *__covar, double t, double dt)\n{\n",
skeleton="\nvoid {%name%}_skelfn (double *__f, double *__x, double *__p, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n"
)
@@ -91,7 +91,7 @@
footer <- list(
rmeasure="\n}\n\n",
- dmeasure="\nlik = (__give_log) ? log(lik) : lik;\n}\n\n",
+ dmeasure="\n}\n\n",
step.fn="\n}\n\n",
skeleton="\n}\n\n"
)
@@ -103,8 +103,8 @@
paramnames <- cleanForC(paramnames)
obsnames <- cleanForC(obsnames)
- modelfile <- paste(name,".c",sep="")
- solib <- paste(name,.Platform$dynlib.ext,sep="")
+ modelfile <- paste0(name,".c")
+ solib <- paste0(name,.Platform$dynlib.ext)
out <- file(description=modelfile,open="w")
@@ -120,10 +120,10 @@
cat(file=out,render(define$var,variable=obsnames[v],ptr='__y',ilist='__obsindex',index=v-1))
}
for (v in seq_along(statenames)) {
- cat(file=out,render(define$var,variable=paste("D",statenames[v],sep=""),ptr='__f',ilist='__stateindex',index=v-1))
+ cat(file=out,render(define$var,variable=paste0("D",statenames[v]),ptr='__f',ilist='__stateindex',index=v-1))
}
cat(file=out,render(define$var.alt,variable="lik",ptr='__lik',index=0))
-
+
## rmeasure function
cat(file=out,render(header$rmeasure,name=name),rmeasure,footer$rmeasure)
@@ -157,14 +157,13 @@
cat(file=out,render(undefine$var,variable=obsnames[v]))
}
for (v in seq_along(statenames)) {
- cat(file=out,render(undefine$var,variable=paste("D",statenames[v],sep="")))
+ cat(file=out,render(undefine$var,variable=paste0("D",statenames[v])))
}
close(out)
- cflags <- paste("PKG_CFLAGS=\"",
+ cflags <- paste0("PKG_CFLAGS=\"",
Sys.getenv("PKG_CFLAGS"),
- " -I",system.file("include",package="pomp"),"\"",
- sep="")
+ " -I",system.file("include",package="pomp"),"\"")
rv <- system2(
command=R.home("bin/R"),
@@ -203,5 +202,5 @@
}
retval[[i]] <- tpl
}
- do.call(function(...)paste(...,sep='.'),retval)
+ do.call(paste0,retval)
}
Modified: pkg/pomp/inst/ChangeLog
===================================================================
--- pkg/pomp/inst/ChangeLog 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/inst/ChangeLog 2012-04-24 17:09:37 UTC (rev 679)
@@ -1,5 +1,28 @@
+2012-04-24 kingaa
+
+ * [r678] R/aaa.R, R/nlf-funcs.R, R/plot-pomp.R, R/pomp.R,
+ R/probe-match.R, R/probe.R, R/spect.R: - use 'paste0' where
+ appropriate. define 'paste0' for use with R <= 2.15.0
+
+2012-04-23 kingaa
+
+ * [r677] DESCRIPTION: - change Depends to 2.14.2
+ * [r676] DESCRIPTION, R/builder.R, data/bbs.rda,
+ data/blowflies.rda, data/dacca.rda, data/euler.sir.rda,
+ data/gillespie.sir.rda, data/gompertz.rda, data/ou2.rda,
+ data/ricker.rda, data/rw2.rda, data/verhulst.rda, inst/NEWS,
+ inst/data-R/blowflies.R, src/blowfly.c, tests/blowflies.R,
+ tests/blowflies.Rout.save, tests/sir.Rout.save: - remove the
+ delay 'tau' from the blowflies example
+ - make the 'pompBuilder' work under 2.14.2
+ - update the data()-loadable examples
+
2012-04-21 kingaa
+ * [r675] inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/intro_to_pomp.pdf: - update vignettes
+ * [r674] DESCRIPTION, inst/ChangeLog, inst/NEWS: - bump version
+ number since R 2.15.0 is now needed
* [r673] inst/NEWS, tests/pfilter.Rout.save: - update unit tests
and NEWS
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/inst/NEWS 2012-04-24 17:09:37 UTC (rev 679)
@@ -1,4 +1,7 @@
NEWS
+0.41-8
+ o A demonstration of 'pompBuilder' has been put into the 'advanced topics' vignette.
+
0.41-7
o A bug in the 'blowflies' example has been fixed.
Thanks to Greg Minshall for discovering it.
Modified: pkg/pomp/inst/TODO
===================================================================
--- pkg/pomp/inst/TODO 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/inst/TODO 2012-04-24 17:09:37 UTC (rev 679)
@@ -1,3 +1,5 @@
+* userdata access from C functions
+
* pompBuilder function for automatic pomp construction
* effects of replacing 'sample.int' with systematic resampling inside 'bsmc'
Modified: pkg/pomp/inst/doc/Makefile
===================================================================
--- pkg/pomp/inst/doc/Makefile 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/inst/doc/Makefile 2012-04-24 17:09:37 UTC (rev 679)
@@ -20,6 +20,7 @@
clean:
$(RM) *.tex *.log *.aux *.blg *.bbl *.out *.Rout *.toc *.lof *.lot
+ $(RM) *.c *.o *.so
$(RM) Rplots.pdf
$(RM) intro_to_pomp-*.pdf advanced_topics_in_pomp-*.pdf
$(RM) intro_to_pomp-*.png advanced_topics_in_pomp-*.png
Modified: pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw 2012-04-24 16:23:25 UTC (rev 678)
+++ pkg/pomp/inst/doc/advanced_topics_in_pomp.Rnw 2012-04-24 17:09:37 UTC (rev 679)
@@ -549,6 +549,143 @@
The \R\ scripts that generated these are included in the \code{data-R} directory of the installed package.
The majority of these use compiled code, which can be found in the package source.
+\section{Pomp Builder}
+
+<<pomp-builder,eval=F>>=
+rmeas <- "
+ double size = 1.0/sigma/sigma;
+ double prob = 1.0/(1.0+rho*cases/size);
+ reports = rnbinom(size,prob);
+"
+
+dmeas <- "
+ double size = 1.0/sigma/sigma;
+ double prob = 1.0/(1.0+rho*cases/size);
+ lik = dnbinom(reports,size,prob,give_log);
+"
+
+stepfn <- "
+ int nrate = 6;
+ int nbasis = 3;
+ int degree = 3;
+ double period = 1.0;
+ double rate[nrate]; // transition rates
+ double trans[nrate]; // transition numbers
+ double dW;
+ double seasonality[nbasis];
+ double beta;
+ int k;
+
+ dW = rgammawn(beta_sd,dt); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
+ periodic_bspline_basis_eval(t,period,degree,nbasis,seasonality);
+ beta = exp(log(beta1)*seasonality[0]+log(beta2)*seasonality[1]+log(beta3)*seasonality[2]);
+
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*I*dW/dt)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the transition numbers
+ trans[0] = rpois(rate[0]*dt); // births are Poisson
+ reulermultinom(2,S,&rate[1],dt,&trans[1]);
+ reulermultinom(2,I,&rate[3],dt,&trans[3]);
+ reulermultinom(1,R,&rate[5],dt,&trans[5]);
+
+ // balance the equations
+ S += trans[0]-trans[1]-trans[2];
+ I += trans[1]-trans[3]-trans[4];
+ R += trans[3]-trans[5];
+ cases += trans[3]; // cases are cumulative recoveries
+ if (beta_sd > 0.0) W += (dW-dt)/beta_sd; // mean = 0, variance = dt
+"
+
+skel <- "
+ int nrate = 6;
+ int nbasis = 3;
+ int degree = 3; // degree of seasonal basis functions
+ double period = 1.0;
+ double rate[nrate]; // transition rates
+ double term[nrate]; // terms in the equations
+ double beta;
+ double seasonality[nbasis];
+
+ // compute transmission rate from seasonality
+ periodic_bspline_basis_eval(t,period,degree,nbasis,seasonality);
+ beta = exp(log(beta1)*seasonality[0]+log(beta2)*seasonality[1]+log(beta3)*seasonality[2]);
+
+ // compute the transition rates
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*I)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the several terms
+ term[0] = rate[0];
+ term[1] = rate[1]*S;
+ term[2] = rate[2]*S;
+ term[3] = rate[3]*I;
+ term[4] = rate[4]*I;
+ term[5] = rate[5]*R;
+
+ // assemble the differential equations
+ DS = term[0]-term[1]-term[2];
+ DI = term[1]-term[3]-term[4];
+ DR = term[3]-term[5];
+ Dcases = term[3]; // accumulate the new I->R transitions
+ DW = 0;
+"
+
+pompBuilder(
+ data=data.frame(
+ reports=NA,
+ time=seq(0,10,by=1/52)
+ ),
+ times="time",
+ t0=-1/52,
+ name="SIR",
+ step.fn.delta.t=1/52/20,
+ paramnames=c(
+ "beta1","beta2","beta3","gamma","mu",
+ "beta.sd","rho","popsize","iota","sigma"
+ ),
+ statenames=c("S","I","R","W","cases"),
+ zeronames="cases",
+ rmeasure=rmeas,
+ dmeasure=dmeas,
+ step.fn=stepfn,
+ skeleton=skel,
+ skeleton.type="vectorfield"
+ ) -> sir
+
+simulate(
+ sir,
+ params=c(gamma=26,mu=0.05,beta.sd=0.1,
+ rho=0.6,sigma=0.1,popsize=1e5,iota=10,
+ beta1=100,beta2=120,beta3=80,
+ S.0=26000,I.0=0,R.0=74000,W.0=0,cases.0=0
+ )
+ ) -> sir
+@
+
+<<pomp-builder-eval,echo=F,eval=T,results=hide>>=
+if (Sys.getenv("POMP_BUILD_VIGNETTES")=="yes") {
+ require(pomp)
+ pompBuilder <- pomp:::pompBuilder
+<<pomp-builder>>
+ simulate(sir) -> x
+ pfilter(sir,Np=500) -> pf
+ dyn.unload("SIR.so")
+ print(logLik(pf))
+}
+@
+
+
+
<<restore-opts,echo=F,results=hide>>=
options(glop)
@
Modified: pkg/pomp/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
More information about the pomp-commits
mailing list