# inst/OldDemo/gillespie.r In CollocInfer: Collocation Inference for Dynamic Systems

```##
## Gillespie simulator for forced SEIR model
##
##    For given parameter values and initial conditions, compute
##    the deterministic SEIR solution and some Gillespie realizations.
##    Plot the results and save to a file (gillespieInc.csv).
##
## Created:  7 Mar 2008 by David Earn
## Changed: 17 Mar 2008
##

compute.gillespie <- TRUE; # FALSE means read gillespie from file

#if (interactive()) {
#} else {
#  # postscript(psfile);
#  pdffile <- "gillespieInc.pdf";
#  pdf(pdffile, width=8.5, height=11);
#}

require(odesolve);  ## for lsoda

## SET PARAMETER VALUES
N <- 5*10^6; ## population size
nu <- 0.02; ## birth rate (1/yr)
mu <- 0.02; ## death rate (1/yr)
R0 <- 17; ## reproductive number
latent.days <- 8;  ## mean latent period (days)
infectious.days <- 5;  ## mean infectious period (days)
alpha <- 0.08; ## seasonal forcing amplitude

## SET INITIAL CONDITIONS
## Chosen to be close to attractor.
s0 <- 0.0576*N;
e0 <- 0.000165*N;
i0 <- 0.0001*N;
if (N>=10^4) {
s0 <- round(s0);
e0 <- round(e0);
i0 <- round(i0);
}

## SET SIMULATION TIME (yrs)
tstart <- 0;
tend <- 5;
tby <- 1/52;
#tend <- 0.05;
#tby <- 0.001;
report.times <- seq(from=tstart, to=tend, by=tby);

## SET NUMBER OF GILLESPIE REALIZATIONS
Nreal <- 1;

## SET DERIVED PARAMETERS
sigma <- 365/latent.days;
gamma <- 365/infectious.days;
beta0 <- R0 * gamma / N; # FIX: this is approximate
Ivisitors <- 0;

## DEFINE PARMS AND IC VECTORS FOR CONVENIENCE
parms <- c(
N=N,
nu=nu,
mu=mu,
R0=R0,
sigma=sigma,
gamma=gamma,
beta0=beta0,
Ivisitors=Ivisitors,
alpha=alpha
);
## we save cumulative incidence as well so we can
## plot incidence when we are done
ic <- c( S=s0, E=e0, I=i0, R=N-s0-e0-i0, cumInc=0 );
cat("Initial conditions:\n");
print(ic);

## SEASONALLY FORCED TRANSMISSION RATE FUNCTION
betafun <- function( t, parms ) {
beta0 <- parms["beta0"];
alpha <- parms["alpha"];
return( beta0 * (1 + alpha*cos(2*pi*t) ))
}

## DERIVATIVE FUNCTION FOR SEIR MODEL
deriv.seir <- function( t, x, parms ) {
with(as.list(x,parms),{
dx <- x; ## dx is same length as x and has same component names
dx["S"] <- nu*N - betafun(t,parms)*S*I - mu*S;
dx["E"] <- betafun(t,parms)*S*I - (sigma+mu)*E;
dx["I"] <- sigma*E - (gamma+mu)*I;
dx["R"] <- gamma*I - mu*R;
## integrate incidence to get cumulative incidence: to get
## incidence we will take adjacent differences of solution
dx["cumInc"] <- betafun(t,parms)*S*I;
return(list(dx))
});
}

## RUN DETERMINISTIC MODEL
lsoda.out <- lsoda( y=ic, times=report.times, func=deriv.seir, parms=parms );

## FUNCTION TO EXTRACT LAST COLUMN OF A MATRIX
lastcol <- function( x ) {
return( x[,ncol(x)] )
}

## PLOT THE DETERMINISTIC SOLUTION
par(mfrow=c(2,1)); #ncol,nrow
t <- lsoda.out[,1];
## extract deterministic cumulative incidence:
detinc <- lastcol(lsoda.out);
## convert to incidence:
detinc <- c(0,diff(detinc));
## plot incidence:
plot(t,detinc,typ="l",col="dark red",lwd=5,xlab="Time (years)",ylab="Incidence",ylim=c(0,1.1*max(detinc)));

#########################
## GILLESPIE ALGORITHM ##
#########################
## Note: This not the most general formulation of the Gillespie
##       algorithm.  More generally we would consider a transition
##       matrix with columns for event names, event rates,  and
##       a column for each state variable.  In each row, the state
##       variable column would contain an integer indicating by
##       how much it increases or decreases as a result of the
##       event in question.  In the SEIR case, these integers
##       would all be -1, 0 or 1 and at most two in any row would
##       be non-zero.

## EVENT RATES
event.rates.seir <- function( t, x, parms ) {
with(as.list(x,parms),{
transmission <- as.vector(betafun(t,parms)*S*(I+Ivisitors));
## as.vector prevents s2e becoming s2e.beta0
rates <- c(snew=nu*N, ## new susceptibles
s2d=mu*S,e2d=mu*E,i2d=mu*I,r2d=mu*R, ## death
s2e=transmission,
e2i=sigma*E, ## becoming infectious
i2r=gamma*I  ## recovery
);
return(rates)
});
}

## GLOBAL VECTORS GIVING COMPARTMENT INCREMENTS AND DECREMENTS
event.names <- c(
"snew", "s2d", "e2d", "i2d", "r2d", "s2e", "e2i", "i2r" );
event.plus <- NULL;
event.plus[event.names] <- c(
"S",    NA,   NA,   NA,   NA,   "E",   "I",   "R" );
event.minus <- NULL;
event.minus[event.names] <- c(
NA,    "S",   "E",   "I",   "R",   "S",   "E",   "I" );
event.counts <- NULL; # vector of counts of each type of event
event.counts[event.names] <- 0;

## TIME TO NEXT GILLESPIE EVENT
time.to.next.event <- function( event.rates ) {
total.rate <- sum(event.rates);
stopifnot(total.rate > 0);
u <- runif(1);
stopifnot(u>0);
return(log(1/u)/total.rate)
##
##
##           return(rexp(1,total.rate))
##
##       which would definitely be preferable if we were
##       dealing with non-exponential distributions that
##       aren't so easy to invert.  However, if the stage
##       durations are not exponentially distributed then
##       the algorithm changes.  In the special case of
##       Erlang distributed stage durations we can stick
##       with the identical algorithm at the expense of
}

## FUNCTION TO RETURN LAST ELEMENT OF AN OBJECT
## Strange that this is not available already in R.
last <- function(x) {
return(x[length(x)])
}

## TYPE OF NEXT GILLESPIE EVENT
## FIX: depends on global vector event.names
select.event.type <- function( event.rates ) {
ne <- length(event.rates);
event.divisions <- c(0,cumsum(event.rates));
event.divisions <- event.divisions/event.divisions[ne+1];
## Note that event.divisions[ne+1] == sum(event.rates)
r <- runif(1);
event.names.index <- 0;
event.names.index <- last(which(event.divisions < r));
stopifnot(event.names.index <= ne);
return(event.names[event.names.index])
}

## UPDATE DISCRETE STATE
## FIX: depends on global vectors event.plus and event.minus
##      (This would be cured by passing a transition matrix.)
update.state <- function(x,event.rates,incidence.count) {
## determine event type
event.name <- select.event.type( event.rates );
## increment and decrement appropriate compartments
if (!is.na(event.plus[event.name])) {
x[event.plus[event.name]] <- x[event.plus[event.name]] + 1;
}
if (!is.na(event.minus[event.name])) {
x[event.minus[event.name]] <- x[event.minus[event.name]] - 1;
}
## keep track of numbers of each event type for later analysis
event.counts[event.name] <<- event.counts[event.name] + 1;
## N.B. <<- is global assign
## INCREMENT INCIDENCE IF THIS WAS AN e2i EVENT:
## FIX: This is an ugly hack... keeping track of a quantity
##      that is not itself a state variable would be easier
##      with a transition matrix formulation.
##      Having added a compartment for incidence in the
##      deterministic model, it would have been better here
##      to increment x["cumInc"] rather than create a new
##      global variable for incidence...
if (event.name == "e2i") {
##cat("incidence.count =", incidence.count, "\n");
incidence.count <<- incidence.count + 1;
}
return(x)
}

## RUN GILLESPIE REALIZATION
gillespie <- function( ic, times, event.rates.func, parms ) {
with(as.list(parms), {
## create matrix in which to return simulation results
nt <- length(times); ## number of time points at which to save
result <- matrix(0,nrow=nt,ncol=1+length(ic));
incidence <- matrix(0,nrow=nt,ncol=2);
## set and save initial state
it <- 1; ## index of time points at which to save
t <- times[it]; ## current time
x <- ic; ## current state
result[it,] <- c(t,x);
## iterate Gillespie algorithm
event.count <- 0;
while (it < nt) {
event.count <- event.count + 1;
if (x["E"]==0 && x["I"]==0) {
cat(sprintf("Extinct at time %g after %d events\n", t, event.count));
break; # Beware: this also stops the birth-death process
}
event.rates <- event.rates.func(t,x,parms);
dt <- time.to.next.event(event.rates);
if (FALSE) { # debugging code
##if (sum(x) > N) {
##if (event.count %% 1000 == 0) {
cat("ireal=", ireal, " it=",it," dt=", dt, " t=", t, " x=", x,
" sum(x)=", sum(x), "\n");
}
t <- t+dt;
if (t > times[it+1]) {
it <- it+1;
result[it,] <- c(times[it],x); # prevalence is saved here
incidence[it,] <- c(times[it],incidence.count);
if (it %% 10 == 0) {
cat("ireal=", ireal, "it=", it, " inc=", incidence[it,2], " SEIR=", result[it,], "\n");
}
incidence.count <<- 0; # restart incidence counter
}
x <- update.state(x,event.rates,incidence.count);
}
## return(result) # return full SEIR at each save point
return(incidence) # return incidence between save points
})
}

#############################
## END GILLESPIE FUNCTIONS ##
#############################

Nvar <- length(ic);
times <- lsoda.out[,1];
all.results <- matrix(0,nrow=length(times),ncol=1+Nvar*(1+Nreal));
all.results[,1] <- times;
all.incidence <- matrix(0,nrow=length(times),ncol=2+Nreal);
all.incidence[,1] <- times;
all.incidence[,2] <- detinc;

outfile <- "gillespieInc.csv";
if (compute.gillespie) {
## COMPUTE SEQUENCE OF GILLESPIE REALIZATIONS
for (ireal in 1:Nreal) {
incidence.count <- 0; # FIX: this is ugly...
gillespie.out <- gillespie(ic=ic,times=report.times,
event.rates.func=event.rates.seir,
parms=parms);
inc <- gillespie.out[,2];
all.incidence[,ireal+2] <- inc;
}
write.csv(all.incidence,outfile,quote=FALSE,row.names=FALSE);
} else {
## READ SAVED GILLESPIE REALIZATIONS FROM FILE
}

## PLOT THE ENSEMBLE MEAN OF THE GILLESPIE REALIZATIONS
t <- all.incidence[,1];
if (Nreal > 1) {
meanreal <- rowMeans(all.incidence[,-c(1,2)]);
} else {
meanreal <- all.incidence[,-c(1,2)];
}
points(t,meanreal,lwd=3);
title(main=sprintf("Determinstic SEIR (red) and mean of %d Gillespie realizations (circles)", Nreal));

## PLOT THE GILLESPIE REALIZATIONS
matplot(t,all.incidence[,-c(1,2)],typ="l",lty=1,
xlab="Time (years)",
ylab="Incidence");
points(t,meanreal,lwd=3);
title(main=sprintf("%d Gillespie realizations (lines) and mean (circles), N = %g", Nreal, N));

## SAVE TIME REQUIRED FOR THIS RUN IN CASE WE WANT TO REDO IT
ptm <- as.vector(proc.time());
## dump proc.time to output file (multiple formats is just
## laziness... should have a function that computes hours,mins,secs
## and outputs in one nice format
if(compute.gillespie){
cat("#      user   system  elapsed\n",
"# seconds:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm, ptm, ptm),
"# minutes:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm/60, ptm/60, ptm/60),
"# hours:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm/60/60, ptm/60/60, ptm/60/60),
file=outfile,append=TRUE);
}
cat("#      user   system  elapsed\n",
"# seconds:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm, ptm, ptm),
"# minutes:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm/60, ptm/60, ptm/60),
"# hours:\n",
sprintf("# %8.2f %8.2f %8.2f\n", ptm/60/60, ptm/60/60, ptm/60/60));

## SAVE EVENT COUNTS FOR ANALYSIS
if (compute.gillespie) {
event.count <- sum(event.counts);
events.per.second <- event.count / ptm;
print(event.counts);
cat(sprintf("# total event.count: %d,  events.per.second: %g\n", event.count, events.per.second));
cat(sprintf("# total event.count: %d,  events.per.second: %g\n", event.count, events.per.second),file=outfile,append=TRUE);
}
```

## Try the CollocInfer package in your browser

Any scripts or data that you put into this service are public.

CollocInfer documentation built on May 2, 2019, 4:03 a.m.