[Stpp-commits] r35 - in pkg: R doc doc/figures man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 21 15:50:52 CEST 2011
Author: gabriele
Date: 2011-07-21 15:50:51 +0200 (Thu, 21 Jul 2011)
New Revision: 35
Added:
pkg/R/STIKhat.r
pkg/R/plotK.R
pkg/R/stan.R
pkg/doc/.dvi
pkg/doc/docV2.aux
pkg/doc/docV2.bbl
pkg/doc/docV2.blg
pkg/doc/docV2.dvi
pkg/doc/docV2.log
pkg/doc/docV2.pdf
pkg/doc/docV2.tex
pkg/doc/figures/fmd.eps
pkg/doc/figures/stik.contour.eps
pkg/doc/figures/stik.persp.eps
pkg/man/STIKhat.Rd
pkg/man/plotK.Rd
pkg/man/stan.Rd
pkg/src/bounds.cmn
pkg/src/stikfunction.f
Removed:
pkg/R/stani.R
pkg/man/stani.Rd
Modified:
pkg/R/covst.r
pkg/R/rinter.r
pkg/R/rlgcp.r
pkg/R/rpp.r
pkg/doc/source.exemples.R
pkg/man/rinfec.Rd
pkg/man/rinter.Rd
pkg/man/rlgcp.Rd
pkg/man/rpcp.Rd
pkg/man/rpp.Rd
pkg/src/stpp.dll
Log:
Added: pkg/R/STIKhat.r
===================================================================
--- pkg/R/STIKhat.r (rev 0)
+++ pkg/R/STIKhat.r 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,110 @@
+STIKhat <- function(xyt, s.region, t.region, dist, times, lambda, correction = TRUE, infectious = TRUE)
+{
+
+# library(splancs)
+ #
+ # This function computes the space-time inhomogeneous K function.
+ #
+ # Arguments:
+ # - xyt: coordinates and times of the point process,
+ # - s.region: set of points defining the domain D,
+ # - t.region: vector giving the lower and upper boundaries of T
+ # - dist: vector of distances h at which K() is computed,
+ # - times: vector of times t at which K() is computed,
+ # - lambda: vector of values of the space-time intensity
+ # function evaluated at the points of the pattern,
+ # - correction: logical value. If True, the Ripley's edge correction
+ # is used.
+ # - infectious: logical value. If True, the STIK function is the one defined
+ # in the case of infectious diseases.
+ #
+ # Value:
+ # - Khat: a ndist x ntimes matrix containing the values of K(h,t),
+ # - dist: the vector of distances used,
+ # - times: the vector of times used.
+ #
+ #
+ # E. Gabriel, september 2005
+ #
+ # last modification: 25.10.2006
+ #
+
+ if (missing(s.region)) s.region <- sbox(xyt[,1:2],xfrac=0.01,yfrac=0.01)
+ if (missing(t.region)) t.region <- range(xyt[,3],na.rm=T)
+
+ pts <- xyt[,1:2]
+ xytimes <- xyt[,3]
+ ptsx <- pts[, 1]
+ ptsy <- pts[, 2]
+ ptst <- xytimes
+ npt <- length(ptsx)
+ ndist <- length(dist)
+ dist <- sort(dist)
+ ntimes <- length(times)
+ times <- sort(times)
+ bsupt <- max(t.region)
+ binft <- min(t.region)
+
+ area <- areapl(s.region)*(bsupt-binft)
+
+ np <- length(s.region[, 1])
+ polyx <- c(s.region[, 1], s.region[1, 1])
+ polyy <- c(s.region[, 2], s.region[1, 2])
+ hkhat <- array(0, dim = c(ndist,ntimes))
+
+ if(missing(lambda))
+ {
+ misl <- 1
+ lambda <- rep(1,npt)
+ }
+ else misl <- 0
+ if (length(lambda)==1) lambda <- rep(lambda,npt)
+ if (correction==T){edg <- 1} else edg <- 0
+ if (infectious==T){infd <- 1} else infd <- 0
+ storage.mode(hkhat) <- "double"
+
+ if (area>10) lambda <- lambda*area
+
+# dyn.load("/home/gabriel/functions/STIKfunction/libF/stikfunction.dll")
+# pathname <- paste(path,"stikfunction.dll",sep="")
+# dyn.load(pathname)
+ nev <- rep(0,ntimes)
+ klist <- .Fortran("stikfunction", as.double(ptsx),
+ as.double(ptsy), as.double(ptst),
+ as.integer(npt), as.double(polyx),
+ as.double(polyy), as.integer(np),
+ as.double(dist), as.integer(ndist),
+ as.double(times), as.integer(ntimes),
+ as.double(bsupt), as.double(binft),
+ as.double(lambda), as.integer(infd),
+ as.integer(edg), (hkhat))
+ khat <- klist[[17]]
+ if (misl==1)
+ {
+ if (area>10)
+ hkhat <- ((area^3)/(npt*(npt-1)))*khat
+ else
+ hkhat <- (area/(npt*(npt-1)))*khat
+ }
+ else
+ {
+ hkhat <- khat
+
+ if (area>10)
+ hkhat <- hkhat*area
+ else
+ hkhat <- hkhat/area
+ }
+ if(dist[1]==0) hkhat[1,]=0
+ if(times[1]==0) hkhat[,1]=0
+
+ Khpp <- matrix(0,ncol=length(times),nrow=length(dist))
+ for(i in 1:length(dist)){Khpp[i,]<-pi*(dist[i]^2)*times}
+
+ if (infectious==T) Ktheo <- Khpp
+ else Ktheo <- 2*Khpp
+
+ lhat <- hkhat-Ktheo
+
+ invisible(return(list(Khat=hkhat,Ktheo=Ktheo,dist=dist,times=times,infectious=infectious)))
+}
Modified: pkg/R/covst.r
===================================================================
--- pkg/R/covst.r 2010-08-31 09:46:16 UTC (rev 34)
+++ pkg/R/covst.r 2011-07-21 13:50:51 UTC (rev 35)
@@ -14,12 +14,12 @@
for (i in 1:length(unique(model)))
{
- if (((isTRUE(separable)) && ((model[i]==models[5]) || (model[i]==models[6]))) || ((!(isTRUE(separable))) && ((model[i]==models[1]) || (model[i]==models[2]) || (model[i]==models[3]) || (model[i]==models[4]) || (model[i]==models[7])))) stop("'stcov' does not match with 'model'")
+ if (((isTRUE(separable)) & ((model[i]==models[5]) | (model[i]==models[6]))) | ((!(isTRUE(separable))) & ((model[i]==models[1]) | (model[i]==models[2]) | (model[i]==models[3]) | (model[i]==models[4]) | (model[i]==models[7])))) stop("'stcov' does not match with 'model'")
}
if (isTRUE(separable))
{
- if ((length(model)!=1) && (length(model)!=2))
+ if ((length(model)!=1) & (length(model)!=2))
stop("for separable covariance functions, 'model' must be of length 1 or 2")
if (length(model)==1)
{
@@ -36,9 +36,9 @@
if (model=="stable")
{
mods <- 2
- if ((param[1] >2) || (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
+ if ((param[1] >2) | (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
modt <- 2
- if ((param[2] >2) || (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
+ if ((param[2] >2) | (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model=="cauchy")
{
@@ -73,12 +73,12 @@
if (model[1]=="stable")
{
mods <- 2
- if ((param[1] >2) || (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
+ if ((param[1] >2) | (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model[2]=="stable")
{
modt <- 2
- if ((param[2] >2) || (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
+ if ((param[2] >2) | (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model[1]=="cauchy")
{
@@ -114,17 +114,17 @@
{
mod <- 5
if (param[6]<2) stop("for Gneiting's covariance function, the sixth parameter must be greater than 2")
- if ((param[3]<=0) || (param[3]>2)) stop("for Gneiting's covariance function, the third parameter must lie in (0,2]")
- if ((param[4]<=0) || (param[4]>1)) stop("for Gneiting's covariance function, the fourth parameter must lie in (0,1]")
- if ((param[5]!=1) && (param[5]!=2) && (param[5]!=3)) stop("for Gneiting's covariance function, the fifth parameter must be 1, 2 or 3")
- if ((param[2]!=1) && (param[2]!=2) && (param[2]!=3)) stop("for Gneiting's covariance function, the second parameter must be 1, 2 or 3")
- if ((param[2]==1) && ((param[1]<0) || (param[1]>2))) stop("for Gneiting's covariance function, if the second parameter equals 1, the first parameter must lie in [0,2]")
- if ((param[2]==2) && (param[1]<=0)) stop("for Gneiting's covariance function, if the second parameter equals 2, the first parameter must be strictly positive")
+ if ((param[3]<=0) | (param[3]>2)) stop("for Gneiting's covariance function, the third parameter must lie in (0,2]")
+ if ((param[4]<=0) | (param[4]>1)) stop("for Gneiting's covariance function, the fourth parameter must lie in (0,1]")
+ if ((param[5]!=1) & (param[5]!=2) & (param[5]!=3)) stop("for Gneiting's covariance function, the fifth parameter must be 1, 2 or 3")
+ if ((param[2]!=1) & (param[2]!=2) & (param[2]!=3)) stop("for Gneiting's covariance function, the second parameter must be 1, 2 or 3")
+ if ((param[2]==1) & ((param[1]<0) | (param[1]>2))) stop("for Gneiting's covariance function, if the second parameter equals 1, the first parameter must lie in [0,2]")
+ if ((param[2]==2) & (param[1]<=0)) stop("for Gneiting's covariance function, if the second parameter equals 2, the first parameter must be strictly positive")
}
if (model=="cesare")
{
mod <- 6
- if (((param[1]>2) || (param[1]<1)) || ((param[2]>2) || (param[2]<1))) stop("for De Cesare's model, the first and second parameters must lie in [1,2]")
+ if (((param[1]>2) | (param[1]<1)) | ((param[2]>2) | (param[2]<1))) stop("for De Cesare's model, the first and second parameters must lie in [1,2]")
if (param[3]<3/2) stop("for De Cesare's model, the third parameter must be greater than 3/2")
}
}
Added: pkg/R/plotK.R
===================================================================
--- pkg/R/plotK.R (rev 0)
+++ pkg/R/plotK.R 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,82 @@
+plotK <- function(K,n=15,L=TRUE,persp=FALSE,legend=TRUE,...)
+{
+ if (isTRUE(L))
+ {
+ k <- K$Khat-K$Ktheo
+ if (isTRUE(K$infectious))
+ titl <- expression(hat(K)[ST] * group("(",list(u,v),")") - pi*u^2*v)
+ else
+ titl <- expression(hat(K)[ST] * group("(",list(u,v),")") - 2*pi*u^2*v)
+ }
+ else
+ {
+ k <- K$Khat
+ titl <- expression(hat(K)[I] * group("(",list(u,v),")") )
+ }
+
+
+ colo <- colorRampPalette(c("red", "white", "blue"))
+ M <- max(abs(range(k)))
+ M <- pretty(c(-M,M),n=n)
+ n <- length(M)
+ COL <- colo(n)
+ if (isTRUE(persp))
+ {
+ mask <- matrix(0,ncol=length(K$times),nrow=length(K$dist))
+ for(i in 1:length(K$dist)){ for(j in 1:length(K$times)){mask[i,j] <- COL[iplace(x=k[i,j],X=M,xinc=diff(range(M))/(n-1))]}}
+ COL <- mask[1:(length(K$dist)-1),1:(length(K$times)-1)]
+
+ if(isTRUE(legend))
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=1,mar=c(0,0,3,0))
+ par(fig=c(0,0.825,0,1))
+ persp(x=K$dist, y=K$times, z=k, xlab="u",ylab="v", zlab="",expand=1, col=COL, ...)
+ title(titl,cex.main=2)
+ par(fig=c(0.825,1,0,1))
+ mini <- iplace(x=min(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ maxi <- iplace(x=max(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
+ }
+ else
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=1)
+ persp(x=K$dist, y=K$times, z=k, xlab="u",ylab="v", zlab="", expand=1, col=COL, ...)
+ title(titl,cex.main=2)
+ }
+ }
+ else
+ {
+ if(isTRUE(legend))
+ {
+ par(cex.lab=1.5,cex.axis=1.5,font=2,lwd=1,plt=c(0,1,0,1),mar=c(0.5,0.5,2.5,0.5),las=1)
+ par(fig=c(0.1,0.825,0.1,1))
+ contour(K$dist, K$times, k, labcex=1.5,levels=M,drawlabels=F,col=colo(n),zlim=range(M),axes=F)
+ box(lwd=2)
+ at <- axTicks(1)
+ axis(1,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(1,at=at[length(at)],labels="u",cex.axis=2)
+ at <- axTicks(2)
+ axis(2,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(2,at=at[length(at)],labels="v",cex.axis=2)
+ title(titl,cex.main=2)
+ par(fig=c(0,1,0.1,1))
+ mini <- iplace(x=min(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ maxi <- iplace(x=max(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
+ }
+ else
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=2,las=1)
+ contour(K$dist, K$times, k, labcex=1.5,levels=M,drawlabels=T,col=colo(n),zlim=range(M),axes=F)
+ box(lwd=2)
+ at <- axTicks(1)
+ axis(1,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(1,at=at[length(at)],labels="u",cex.axis=2)
+ at <- axTicks(2)
+ axis(2,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(2,at=at[length(at)],labels="v",cex.axis=2)
+ title(titl,cex.main=2)
+ }
+ }
+}
+
Modified: pkg/R/rinter.r
===================================================================
--- pkg/R/rinter.r 2010-08-31 09:46:16 UTC (rev 34)
+++ pkg/R/rinter.r 2011-07-21 13:50:51 UTC (rev 35)
@@ -396,10 +396,10 @@
if (!(is.function(hs)))
{
- if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && (npoints * pi * deltas^2/4 > areapl(s.region)))
+ if ((inhibition==TRUE) & (thetas==0) & (hs=="step") & (npoints * pi * deltas^2/4 > areapl(s.region)))
stop(paste("s.region is too small to fit", npoints, "points", "at minimum distance", deltas))
- if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && ((max(t.region)-min(t.region))/deltat<npoints))
+ if ((inhibition==TRUE) & (thetas==0) & (hs=="step") & ((max(t.region)-min(t.region))/deltat<npoints))
stop(paste("t.region is too small to fit", npoints, "points", "at minimum time interval", deltat))
}
Modified: pkg/R/rlgcp.r
===================================================================
--- pkg/R/rlgcp.r 2010-08-31 09:46:16 UTC (rev 34)
+++ pkg/R/rlgcp.r 2011-07-21 13:50:51 UTC (rev 35)
@@ -177,7 +177,7 @@
npts <- round(lambdamax/(s.area*t.area),0)
if (npts==0) stop("there is no data to thin")
- if ((replace==FALSE) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
+ if ((replace==FALSE) & (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==TRUE)
{
Modified: pkg/R/rpp.r
===================================================================
--- pkg/R/rpp.r 2010-08-31 09:46:16 UTC (rev 34)
+++ pkg/R/rpp.r 2011-07-21 13:50:51 UTC (rev 35)
@@ -229,7 +229,7 @@
x <- xy[,1]
y <- xy[,2]
- if ((replace==F) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
+ if ((replace==F) & (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==T)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
@@ -327,7 +327,7 @@
npts <- npoints
if (npts==0) stop("there is no data to thin")
- if ((replace==F) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
+ if ((replace==F) & (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==T)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
Added: pkg/R/stan.R
===================================================================
--- pkg/R/stan.R (rev 0)
+++ pkg/R/stan.R 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,220 @@
+
+.listmerge = function (x, y, ...)
+{
+### taken from RCurl
+ if (length(x) == 0)
+ return(y)
+ if (length(y) == 0)
+ return(x)
+ i = match(names(y), names(x))
+ i = is.na(i)
+ if (any(i))
+ x[names(y)[which(i)]] = y[which(i)]
+ x
+}
+
+
+.stan3d.redraw <- function(o) {
+ ## switch off redraws
+ par3d(skipRedraw=TRUE)
+ np=dim(o$xyt)[1]
+
+ ## compute new states
+ tin = rep(1,np)
+ tin[o$xyt[,3]>(o$t-o$width)]=2
+ tin[o$xyt[,3]>o$t]=3
+
+ ## which points have changed state since last time?
+ changed = tin != o$xyt[,5]
+
+ ## remove any that changed
+ if(any(changed)){
+ rgl.pop(id=o$xyt[changed,4])
+ }
+
+ ## now add them back in their correct state:
+ for(i in (1:np)[changed]){
+
+### should be as simple as this:
+### material3d(o$states[[tin[i]]])
+### but setting alpha is causing problems. Bug reported. Hence:
+
+ sphereList = list(x=o$xyt[i,1],y=o$xyt[i,2],z=o$xyt[i,3],radius=o$states[[tin[i]]]$radius)
+ materialList = o$states[[tin[i]]]
+ pList = .listmerge(sphereList,materialList)
+ o$xyt[i,4]=do.call(spheres3d,pList)
+
+ o$xyt[i,5]=tin[i]
+ }
+ ## start drawing again:
+ par3d(skipRedraw=FALSE)
+
+ ## update the slider positions
+ .store(o)
+
+ ## and return the modified object:
+ return(o)
+}
+
+.rp.stan3d <- function(xyt,tlim,twid,states) {
+ t=tlim[1];width=twid
+ e=new.env()
+ stan.panel <- rp.control(title="space-time animation",
+ xyt=xyt, t=tlim[1], width=twid,
+ states=states,
+ e=e
+ )
+ rp.slider(stan.panel, t, title = "time", from=tlim[1], to=tlim[2], action = .stan3d.redraw,showvalue=TRUE)
+ rp.slider(stan.panel, width, title = "window", from=0, to=diff(tlim), action = .stan3d.redraw,showvalue=TRUE)
+ rp.button(stan.panel,action=function(p){par3d(FOV=0,userMatrix = rotationMatrix(0, 1,0,0));return(p)},title="reset axes")
+ rp.button(stan.panel,action=.store,title="quit",quitbutton=TRUE)
+ rp.do(stan.panel, .stan3d.redraw)
+ rp.block(stan.panel)
+ return(e)
+}
+
+.store = function(panel){
+ assign("t",panel$t,env=panel$e)
+ assign("width",panel$width,env=panel$e)
+ return(panel)
+}
+
+stan <- function(xyt,tlim=range(xyt[,3],na.rm=TRUE),twid=diff(tlim)/20,persist=FALSE,states,bgpoly,bgframe=TRUE,bgimage,bgcol=gray(seq(0,1,len=12)),axes=TRUE){
+ require(rgl)
+ require(rpanel)
+
+X=(xyt[,1]-min(xyt[,1]))/(diff(range(xyt[,1])))
+Y=(xyt[,2]-min(xyt[,2]))/(diff(range(xyt[,2])))
+Z=(xyt[,3]-min(xyt[,3]))/(diff(range(xyt[,3])))
+
+tlim=(tlim-min(xyt[,3]))/(diff(range(xyt[,3])))
+
+if(!missing(bgpoly)) bgpoly=cbind((bgpoly[,1]-min(xyt[,1]))/(diff(range(xyt[,1]))),(bgpoly[,2]-min(xyt[,2]))/(diff(range(xyt[,2]))))
+
+
+XYT=as.3dpoints(X,Y,Z)
+
+ if(missing(states)){
+ ## default colouring scheme:
+ states=list(
+ past=list(col="blue",radius=1/80,alpha=0.5,lit=FALSE),
+ present=list(col="red",radius=1/30,alpha=0.5,lit=FALSE),
+ ## still-to-come points are invisible (alpha=0)
+ future=list(col="yellow",alpha=0.0,radius=1/80,lit=FALSE)
+ )
+ if(persist){
+ states$past=states$present
+ }
+ }
+
+ maxRadius = max(states$past$radius,states$present$radius,states$future$radius)
+ xr = range(XYT[,1],na.rm=TRUE)
+ yr = range(XYT[,2],na.rm=TRUE)
+ tr = range(XYT[,3],na.rm=TRUE)
+ diag=sqrt(diff(xr)^2+diff(yr)^2+diff(tr)^2)
+
+ states$past$radius=states$past$radius*diag
+ states$present$radius=states$present$radius*diag
+ states$future$radius=states$future$radius*diag
+
+ .setPlot(xr[1],xr[2],yr[1],yr[2],tr[1],tr[2],maxRadius,xyt=xyt,axes=axes)
+
+ if(!missing(bgpoly)){
+ poly=rbind(bgpoly,bgpoly[1,])
+ poly=cbind(poly,min(tr))
+ lines3d(poly,size=2.0)
+ if(bgframe){
+ ci=chull(bgpoly)
+ nci=length(ci)
+ cpoints=bgpoly[ci,]
+ cpoints2 = cpoints[rep(1:nci,rep(2,nci)),]
+ cpoints2 = cbind(cpoints2,c(min(tr),max(tr)))
+ segments3d(cpoints2)
+ poly=cbind(poly[,1:2],max(tr))
+ lines3d(poly,size=2.0)
+ }
+ }
+
+ if(!missing(bgimage)){
+ .setBG(bgimage,min(tr),col=bgcol)
+ }
+
+ XYT=data.frame(XYT[,1:3])
+ XYT$id=NA
+ ## initially all points will need redrawing:
+ XYT$state=-1
+
+
+ ## these points will get redrawn immediately... probably a better way to do this:
+ cat("Setting up...")
+ npts = dim(XYT)[1]
+ if(npts>=100){
+ tenths = as.integer(seq(1,npts,len=10))
+ }
+
+ for(i in 1:(dim(XYT)[1])){
+ if(npts>=100){
+ tn = (1:10)[tenths == i]
+ if(length(tn)>0){
+ cat(paste(11-tn,"...",sep=""))
+ }
+ }
+ XYT[i,4]=points3d(XYT[i,1,drop=FALSE],XYT[i,2,drop=FALSE],XYT[i,3,drop=FALSE],alpha=0.0)
+ }
+ cat("...done\n")
+ env = .rp.stan3d(XYT,tlim,twid,states)
+ ret=list()
+ for(n in ls(env)){
+ ret[[n]]=get(n,env=env)
+ }
+ return(ret)
+}
+
+.ranger <- function(x,margin=0.2){
+ lim=range(x,na.rm=TRUE)
+ lim=lim+c(-margin,margin)*diff(lim)
+ return(lim)
+}
+
+.setPlot=function(xmin,xmax,ymin,ymax,tmin,tmax,radius=1/20,xyt,axes){
+ require(rgl)
+ diag=sqrt((xmax-xmin)^2+(ymax-ymin)^2+(tmax-tmin)^2)
+ xr=c(xmax,xmin)+c(radius*(xmax-xmin),-radius*(xmax-xmin))*2
+ yr=c(ymax,ymin)+c(radius*(ymax-ymin),-radius*(ymax-ymin))*2
+ tr=c(tmax,tmin)+c(radius*(tmax-tmin),-radius*(tmax-tmin))*2
+
+ plot3d(xr,yr,tr,type="n",col="red",box=TRUE,axes=FALSE,xlab="x",ylab="y",zlab="t")
+
+ if(axes==TRUE)
+ {
+ xr2=min(xyt[,1])+range(xr)*(diff(range(xyt[,1])))
+ yr2=min(xyt[,2])+range(yr)*(diff(range(xyt[,2])))
+ tr2=min(xyt[,3])+range(tr)*(diff(range(xyt[,3])))
+ xl=round(seq(min(xr2),max(xr2),length=5))
+ yl=round(seq(min(yr2),max(yr2),length=5))
+ tl=round(seq(min(tr2),max(tr2),length=5))
+
+ axis3d('x-',at=seq(0,1,length=5),labels=xl,nticks=5)
+ axis3d('y-',at=seq(0,1,length=5),labels=yl,nticks=5)
+ axis3d('z-',at=seq(0,1,length=5),labels=tl,nticks=5)
+ par3d(FOV=0)
+ AR=(xmax-xmin)/(ymax-ymin)
+ aspect3d(AR,1,1)
+ par3d(userMatrix = rotationMatrix(0, 1,0,0))
+ }
+ else
+ {
+ axis3d('x-',labels="")
+ axis3d('y-',labels="")
+ axis3d('z-',labels="")
+ par3d(FOV=0)
+ AR=(xmax-xmin)/(ymax-ymin)
+ aspect3d(AR,1,1)
+ par3d(userMatrix = rotationMatrix(0, 1,0,0))
+ }
+}
+
+.setBG=function(xyz,zplane,col=heat.colors(12)){
+ cols = col[as.integer(cut(xyz$z,length(col)))]
+ surface3d(xyz$x,xyz$y,rep(zplane,prod(dim(xyz$z))),col=cols,lit=FALSE)
+}
Deleted: pkg/R/stani.R
===================================================================
--- pkg/R/stani.R 2010-08-31 09:46:16 UTC (rev 34)
+++ pkg/R/stani.R 2011-07-21 13:50:51 UTC (rev 35)
@@ -1,189 +0,0 @@
-
-.listmerge = function (x, y, ...)
-{
-### taken from RCurl
- if (length(x) == 0)
- return(y)
- if (length(y) == 0)
- return(x)
- i = match(names(y), names(x))
- i = is.na(i)
- if (any(i))
- x[names(y)[which(i)]] = y[which(i)]
- x
-}
-
-
-.stan3d.redraw <- function(o) {
- ## switch off redraws
- par3d(skipRedraw=TRUE)
- np=dim(o$xyt)[1]
-
- ## compute new states
- tin = rep(1,np)
- tin[o$xyt[,3]>(o$t-o$width)]=2
- tin[o$xyt[,3]>o$t]=3
-
- ## which points have changed state since last time?
- changed = tin != o$xyt[,5]
-
- ## remove any that changed
- if(any(changed)){
- rgl.pop(id=o$xyt[changed,4])
- }
-
- ## now add them back in their correct state:
- for(i in (1:np)[changed]){
-
-### should be as simple as this:
-### material3d(o$states[[tin[i]]])
-### but setting alpha is causing problems. Bug reported. Hence:
-
- sphereList = list(x=o$xyt[i,1],y=o$xyt[i,2],z=o$xyt[i,3],radius=o$states[[tin[i]]]$radius)
- materialList = o$states[[tin[i]]]
- pList = .listmerge(sphereList,materialList)
- o$xyt[i,4]=do.call(spheres3d,pList)
-
- o$xyt[i,5]=tin[i]
- }
- ## start drawing again:
- par3d(skipRedraw=FALSE)
-
- ## update the slider positions
- .store(o)
-
- ## and return the modified object:
- return(o)
-}
-
-.rp.stan3d <- function(xyt,tlim,twid,states) {
- t=tlim[1];width=twid
- e=new.env()
- stan.panel <- rp.control(title="space-time animation",
- xyt=xyt, t=tlim[1], width=twid,
- states=states,
- e=e
- )
- rp.slider(stan.panel, t, title = "time", from=tlim[1], to=tlim[2], action = .stan3d.redraw,showvalue=TRUE)
- rp.slider(stan.panel, width, title = "window", from=0, to=diff(tlim), action = .stan3d.redraw,showvalue=TRUE)
- rp.button(stan.panel,action=function(p){par3d(FOV=0,userMatrix = rotationMatrix(0, 1,0,0));return(p)},title="reset axes")
- rp.button(stan.panel,action=.store,title="quit",quitbutton=TRUE)
- rp.do(stan.panel, .stan3d.redraw)
- rp.block(stan.panel)
- return(e)
-}
-
-.store = function(panel){
- assign("t",panel$t,env=panel$e)
- assign("width",panel$width,env=panel$e)
- return(panel)
-}
-
-stani <- function(xyt,tlim=range(xyt[,3],na.rm=TRUE),twid=diff(tlim)/20,persist=FALSE,states,bgpoly,bgframe=TRUE,bgimage,bgcol=gray(seq(0,1,len=12))){
- require(rgl)
- require(rpanel)
- if(missing(states)){
- ## default colouring scheme:
- states=list(
- past=list(col="blue",radius=1/80,alpha=0.5,lit=FALSE),
- present=list(col="red",radius=1/30,alpha=0.5,lit=FALSE),
- ## still-to-come points are invisible (alpha=0)
- future=list(col="yellow",alpha=0.0,radius=1/80,lit=FALSE)
- )
- if(persist){
- states$past=states$present
- }
- }
-
- maxRadius = max(states$past$radius,states$present$radius,states$future$radius)
- xr = range(xyt[,1],na.rm=TRUE)
- yr = range(xyt[,2],na.rm=TRUE)
- tr = range(xyt[,3],na.rm=TRUE)
- diag=sqrt(diff(xr)^2+diff(yr)^2+diff(tr)^2)
-
- states$past$radius=states$past$radius*diag
- states$present$radius=states$present$radius*diag
- states$future$radius=states$future$radius*diag
-
- .setPlot(xr[1],xr[2],yr[1],yr[2],tr[1],tr[2],maxRadius)
-
- if(!missing(bgpoly)){
- poly=rbind(bgpoly,bgpoly[1,])
- poly=cbind(poly,min(tr))
- lines3d(poly,size=2.0)
- if(bgframe){
- ci=chull(bgpoly)
- nci=length(ci)
- cpoints=bgpoly[ci,]
- cpoints2 = cpoints[rep(1:nci,rep(2,nci)),]
- cpoints2 = cbind(cpoints2,c(min(tr),max(tr)))
- segments3d(cpoints2)
- poly=cbind(poly[,1:2],max(tr))
- lines3d(poly,size=2.0)
- }
- }
-
- if(!missing(bgimage)){
- .setBG(bgimage,min(tr),col=bgcol)
- }
-
- xyt=data.frame(xyt[,1:3])
- xyt$id=NA
- ## initially all points will need redrawing:
- xyt$state=-1
-
-
- ## these points will get redrawn immediately... probably a better way to do this:
- cat("Setting up...")
- npts = dim(xyt)[1]
- if(npts>=100){
- tenths = as.integer(seq(1,npts,len=10))
- }
-
- for(i in 1:(dim(xyt)[1])){
- if(npts>=100){
- tn = (1:10)[tenths == i]
- if(length(tn)>0){
- cat(paste(11-tn,"...",sep=""))
- }
- }
- xyt[i,4]=points3d(xyt[i,1,drop=FALSE],xyt[i,2,drop=FALSE],xyt[i,3,drop=FALSE],alpha=0.0)
- }
- cat("...done\n")
- env = .rp.stan3d(xyt,tlim,twid,states)
- ret=list()
- for(n in ls(env)){
- ret[[n]]=get(n,env=env)
- }
- return(ret)
-}
-
-.ranger <- function(x,margin=0.2){
- lim=range(x,na.rm=TRUE)
- lim=lim+c(-margin,margin)*diff(lim)
- return(lim)
-}
-
-.setPlot=function(xmin,xmax,ymin,ymax,tmin,tmax,radius=1/20){
- require(rgl)
- diag=sqrt((xmax-xmin)^2+(ymax-ymin)^2+(tmax-tmin)^2)
- xr=c(xmax,xmin)+c(radius*(xmax-xmin),-radius*(xmax-xmin))*2
- yr=c(ymax,ymin)+c(radius*(ymax-ymin),-radius*(ymax-ymin))*2
- tr=c(tmax,tmin)+c(radius*(tmax-tmin),-radius*(tmax-tmin))*2
-
-
- plot3d(xr,yr,tr,type="n",col="red",box=TRUE,axes=FALSE,xlab="x",ylab="y",zlab="t")
- axis3d('x-',tick=FALSE)
- axis3d('y-',tick=FALSE)
- axis3d('z-')
- par3d(FOV=0)
- AR=(xmax-xmin)/(ymax-ymin)
- aspect3d(AR,1,1)
- par3d(userMatrix = rotationMatrix(0, 1,0,0))
-
-}
-
-.setBG=function(xyz,zplane,col=heat.colors(12)){
- cols = col[as.integer(cut(xyz$z,length(col)))]
- surface3d(xyz$x,xyz$y,rep(zplane,prod(dim(xyz$z))),col=cols,lit=FALSE)
-}
Added: pkg/doc/.dvi
===================================================================
(Binary files differ)
Property changes on: pkg/doc/.dvi
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/doc/docV2.aux
===================================================================
--- pkg/doc/docV2.aux (rev 0)
+++ pkg/doc/docV2.aux 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,57 @@
+\relax
+\@writefile{toc}{\contentsline {section}{\numberline {1}Introduction}{1}}
+\@writefile{toc}{\contentsline {section}{\numberline {2}Spatio-Temporal Point Processes}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}First-order and second-order properties}{2}}
+\newlabel{subsec:intensity}{{2.1}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.2}Stationarity}{3}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.3}Separability}{3}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.4}Static and dynamic plotting of spatio-temporal point process data}{3}}
+\@writefile{lof}{\contentsline {figure}{\numberline {1}{\ignorespaces Static two-panel plot of data from the 2001 UK FMD epidemic in the county of Cumbria.}}{4}}
+\newlabel{fig:fmd1}{{1}{4}}
+\@writefile{lof}{\contentsline {figure}{\numberline {2}{\ignorespaces Static plot of data from the 2001 UK FMD epidemic. Time is treated as a quantitative mark; light grey/small dots correspond to the oldest events and dark grey/large dots correspond to the most recent events.}}{4}}
+\newlabel{fig:fmd2}{{2}{4}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.5}Space-time inhomogeneous $K$-function}{5}}
+\@writefile{lof}{\contentsline {figure}{\numberline {3}{\ignorespaces Contour plot (left) and perspective plot (right) of the STIK function estimated from FMD data.}}{6}}
+\newlabel{fig:stik}{{3}{6}}
+\@writefile{toc}{\contentsline {section}{\numberline {3}Models}{7}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}Poisson process}{7}}
+\newlabel{sec:Poisson}{{3.1}{7}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.2}Poisson cluster process}{9}}
+\newlabel{sec:pcp}{{3.2}{9}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.3}Interaction processes}{10}}
+\newlabel{sec:inhib}{{3.3}{10}}
+\@writefile{lof}{\contentsline {figure}{\numberline {4}{\ignorespaces Inhibition functions implemented in {\tt rinter}; {\it left}: step, {\it right}: gaussian.}}{12}}
+\newlabel{fig:hinhib}{{4}{12}}
+\@writefile{lof}{\contentsline {figure}{\numberline {5}{\ignorespaces Contagious functions implemented in {\tt rinter}; {\it left}: step, {\it right}: gaussian.}}{13}}
+\newlabel{fig:hcont}{{5}{13}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.4}Infectious processes}{14}}
+\@writefile{lof}{\contentsline {figure}{\numberline {6}{\ignorespaces Illustration of $\mu (t)$ for $h$ as a step function (left) or a gaussian function (right).}}{16}}
+\newlabel{fig:cuminfec}{{6}{16}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.5}Log-Gaussian Cox processes}{17}}
+\newlabel{eq:sepcov}{{1}{18}}
+\newlabel{eq:GneitingCov}{{2}{18}}
+\@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces {\em Some classes of isotropic covariance functions.}}}{19}}
+\newlabel{tab:modcov}{{1}{19}}
+\@writefile{lot}{\contentsline {table}{\numberline {2}{\ignorespaces {\em Range of the parameters defining Equation (2\hbox {}).}}}{19}}
+\newlabel{tab:nsst}{{2}{19}}
+\bibcite{baddeley2000}{1}
+\bibcite{baddeley2005}{2}
+\bibcite{becker}{3}
+\bibcite{berman1989}{4}
+\bibcite{chan1997}{5}
+\bibcite{cressie1991}{6}
+\bibcite{cox1955}{7}
+\bibcite{cox1980}{8}
+\bibcite{daley2003}{9}
+\bibcite{diggle1985}{10}
+\bibcite{diggle2003}{11}
+\bibcite{diggle2006}{12}
+\bibcite{diggle2010}{13}
+\bibcite{gabriel2009}{14}
+\bibcite{moller2003}{15}
+\bibcite{moller2011}{16}
+\bibcite{neyman1958}{17}
+\bibcite{rowlingson1993}{18}
+\bibcite{silverman1986}{19}
+\bibcite{tanemura1979}{20}
+\bibcite{zhuang2002}{21}
Added: pkg/doc/docV2.blg
===================================================================
--- pkg/doc/docV2.blg (rev 0)
+++ pkg/doc/docV2.blg 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,5 @@
+This is BibTeX, Version 0.99cThe top-level auxiliary file: docV2.aux
+I found no \citation commands---while reading file docV2.aux
+I found no \bibdata command---while reading file docV2.aux
+I found no \bibstyle command---while reading file docV2.aux
+(There were 3 error messages)
Added: pkg/doc/docV2.dvi
===================================================================
(Binary files differ)
Property changes on: pkg/doc/docV2.dvi
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/doc/docV2.log
===================================================================
--- pkg/doc/docV2.log (rev 0)
+++ pkg/doc/docV2.log 2011-07-21 13:50:51 UTC (rev 35)
@@ -0,0 +1,290 @@
+This is pdfTeX, Version 3.1415926-1.40.8-alpha-20080323 (MiKTeX 2.7) (preloaded format=latex 2008.7.6) 21 JUL 2011 15:48
+entering extended mode
+**C:/Documents*and*Settings/edith/Mes*documents/stppWorking/pkg/doc/docV2.tex
+("C:/Documents and Settings/edith/Mes documents/stppWorking/pkg/doc/docV2.tex"
+LaTeX2e <2005/12/01>
+Babel <v3.8j> and hyphenation patterns for english, dumylang, nohyphenation, ge
+rman, ngerman, french, loaded.
+("C:\Program Files\MiKTeX 2.7\tex\latex\base\article.cls"
+Document Class: article 2005/09/16 v1.4f Standard LaTeX document class
+("C:\Program Files\MiKTeX 2.7\tex\latex\base\size11.clo"
+File: size11.clo 2005/09/16 v1.4f Standard LaTeX file (size option)
+)
+\c at part=\count79
+\c at section=\count80
+\c at subsection=\count81
+\c at subsubsection=\count82
+\c at paragraph=\count83
+\c at subparagraph=\count84
+\c at figure=\count85
+\c at table=\count86
+\abovecaptionskip=\skip41
+\belowcaptionskip=\skip42
+\bibindent=\dimen102
+)
+("C:\Program Files\MiKTeX 2.7\tex\latex\fancyhdr\fancyhdr.sty"
+\fancy at headwidth=\skip43
+\f at ncyO@elh=\skip44
+\f at ncyO@erh=\skip45
+\f at ncyO@olh=\skip46
+\f at ncyO@orh=\skip47
+\f at ncyO@elf=\skip48
+\f at ncyO@erf=\skip49
+\f at ncyO@olf=\skip50
+\f at ncyO@orf=\skip51
+)
+("C:\Program Files\MiKTeX 2.7\tex\latex\tools\verbatim.sty"
+Package: verbatim 2003/08/22 v1.5q LaTeX2e package for verbatim enhancements
+\every at verbatim=\toks14
+\verbatim at line=\toks15
+\verbatim at in@stream=\read1
+)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/stpp -r 35
More information about the Stpp-commits
mailing list