[Stpp-commits] r31 - in pkg: . R doc doc/figures man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 13 14:49:52 CET 2010
Author: gabriele
Date: 2010-02-13 14:49:51 +0100 (Sat, 13 Feb 2010)
New Revision: 31
Added:
pkg/doc/
pkg/doc/doc.aux
pkg/doc/doc.dvi
pkg/doc/doc.log
pkg/doc/doc.pdf
pkg/doc/doc.tex
pkg/doc/figures/
pkg/doc/figures/contgauss.eps
pkg/doc/figures/contstep.eps
pkg/doc/figures/cumulative.gaussian.ps
pkg/doc/figures/cumulative.step.ps
pkg/doc/figures/fmdt.eps
pkg/doc/figures/fmdxy.eps
pkg/doc/figures/infection.gaussian.ps
pkg/doc/figures/infection.step.ps
pkg/doc/figures/inhibgauss.eps
pkg/doc/figures/inhibstep.eps
pkg/doc/source.exemples.R
pkg/man/as.3dpoint.Rd
pkg/src/pkg_res.rc
pkg/src/stpp.dll
Removed:
pkg/R/gauss3D.r
pkg/R/iplace.r
pkg/R/make.grid.r
pkg/R/pcp.larger.region.r
pkg/R/rhpp.r
pkg/R/ripp.r
pkg/R/set.cov.r
pkg/R/spatial.inhibition.r
pkg/R/temporal.inhibition.r
pkg/man/as.3dpoints.RD
Modified:
pkg/R/covst.r
pkg/R/is.3dpoints.R
pkg/R/plot.stpp.R
pkg/R/rinter.r
pkg/R/rlgcp.r
pkg/R/rpcp.r
pkg/R/rpp.r
pkg/R/sim.stpp.r
pkg/man/animation.Rd
pkg/man/is.3dpoints.RD
pkg/man/plot.stpp.RD
pkg/man/rinfec.Rd
pkg/man/rinter.Rd
pkg/man/rlgcp.Rd
pkg/man/rpp.Rd
pkg/man/stani.Rd
Log:
Modified: pkg/R/covst.r
===================================================================
--- pkg/R/covst.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/covst.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,3 +1,123 @@
+set.cov <- function(separable,model,param,sigma2)
+ {
+ mods <- 0
+ modt <- 0
+ mod <- 0
+
+ models <- c("exponential","cauchy","stable","wave","gneiting","cesare","matern","none")
+
+ for(i in 1:length(model))
+ {
+ M <- which(models==model[i])
+ if (length(M)==0) stop("the model is not implemented")
+ }
+
+ 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))
+ {
+ if ((length(model)!=1) && (length(model)!=2))
+ stop("for separable covariance functions, 'model' must be of length 1 or 2")
+ if (length(model)==1)
+ {
+ if (model=="none")
+ {
+ mods <- 0
+ modt <- 0
+ }
+ if (model=="exponential")
+ {
+ mods <- 1
+ modt <- 1
+ }
+ if (model=="stable")
+ {
+ mods <- 2
+ if ((param[1] >2) || (param[1]<0)) stop("parameter of stable model must lie in [0,2]")
+ modt <- 2
+ if ((param[2] >2) || (param[2]<0)) stop("parameter of stable model must lie in [0,2]")
+ }
+ if (model=="cauchy")
+ {
+ mods <- 3
+ if (param[1]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ modt <- 3
+ if (param[2]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ }
+ if (model=="wave")
+ {
+ mods <- 4
+ modt <- 4
+ }
+ if (model=="matern") stop("You must specify another model for the temporal covariance")
+ }
+ if (length(model)==2)
+ {
+ if (model[1]=="none")
+ mods <- 0
+ if (model[2]=="none")
+ modt <- 0
+ if (model[1]=="exponential")
+ mods <- 1
+ if (model[2]=="exponential")
+ modt <- 1
+ if (model[1]=="stable")
+ {
+ mods <- 2
+ if ((param[1] >2) || (param[1]<0)) stop("parameter of stable model must lie in [0,2]")
+ }
+ if (model[2]=="stable")
+ {
+ modt <- 2
+ if ((param[2] >2) || (param[2]<0)) stop("parameter of stable model must lie in [0,2]")
+ }
+ if (model[1]=="cauchy")
+ {
+ mods <- 3
+ if (param[1]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ }
+ if (model[2]=="cauchy")
+ {
+ modt <- 3
+ if (param[2]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ }
+ if (model[1]=="wave")
+ mods <- 4
+ if (model[2]=="wave")
+ modt <- 4
+ if (model[1]=="matern")
+ mods <- 7
+ }
+ }
+ if (!(isTRUE(separable)))
+ {
+ if (length(model)!=1)
+ stop("for non-separable covariance functions, 'model' must be of length 1")
+ if (model=="gneiting")
+ {
+ 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 (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[3]<3/2) stop("for De Cesare's model, the third parameter must be greater than 3/2")
+ }
+ }
+
+ return(model=c(mods,modt,mod))
+ }
+
covst <- function(dist,times,separable=TRUE,model,param=c(1,1,1,1,1,2),sigma2=1,scale=c(1,1),plot=TRUE,nlevels=10)
{
@@ -4,7 +124,7 @@
nt <- length(times)
np <- length(dist)
- model <- set.cov(separable,model,param,var.grf)
+ model <- set.cov(separable,model,param,sigma2)
gs <- array(0, dim = c(np,nt))
storage.mode(gs) <- "double"
Deleted: pkg/R/gauss3D.r
===================================================================
--- pkg/R/gauss3D.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/gauss3D.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,78 +0,0 @@
-gauss3D <- function(nx=100,ny=100,nt=100,xlim,ylim,tlim,separable=TRUE,model="exponential",param=c(1,1,1,1,1,2),scale=c(1,1),var.grf=1,mean.grf=0,exact=TRUE)
-{
-
- N <- c(nx,ny,nt)
-
- mod <- set.cov(separable,model,param,var.grf)
-
- g <- floor(log(2*(N-1))/log(2))+1
- M <- 2^g
-
- count <- 1
- changG <- TRUE
- while(changG==TRUE)
- {
- L <- rep(-9999,M[1]*M[2]*M[3])
-
- storage.mode(L) <- "double"
-
- res <- .Fortran("circ",
- (L),
- as.integer(M),
- as.integer(N),
- as.double(xlim),
- as.double(ylim),
- as.double(tlim),
- as.integer(mod),
- as.double(param),
- as.double(var.grf),
- as.double(scale))[[1]]
-
- L <- array(res,dim=M)
- FTL <- fft(L)
-
- if (isTRUE(exact))
- {
- if (min(Re(FTL))<0)
- {
- g <- g+1
- M <- 2^g
- changG <- TRUE
- count <- count+1
- }
- else
- changG <- FALSE
- }
- else
- {
- FTL[Re(FTL)<0] <- 0
- changG <- FALSE
- }
-
- }
- print(count)
-
- X <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
- X <- fft(X,inverse=TRUE)
- A <- sqrt(FTL)*X
- G <- (Re(fft(A,inverse=FALSE))/(M[1]*M[2]*M[3]))[1:N[1],1:N[2],1:N[3]]
- G <- G+mean.grf
-
-## Remarks
-# --------
-#
-# The right way is:
-# X1 <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
-# X2 <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
-# X <- fft(X1+1i*X2,inverse=TRUE)
-# A <- sqrt(FTL)*X
-# G <- (Re(fft(A,inverse=FALSE))/(M[1]*M[2]*M[3]))[1:N[1],1:N[2],1:N[3]]
-# but it doesn't change the results and it is slightly faster.
-#
-# Taking the first elements of FTL and then computing G (changing M by N)
-# is faster than taking the first elements of G, but it does not provide
-# the same results
-
- return(G)
-}
-
Deleted: pkg/R/iplace.r
===================================================================
--- pkg/R/iplace.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/iplace.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,14 +0,0 @@
-iplace <- function(X,x,xinc)
- {
- n <- length(X)
- i <- 0
- repeat
- {
- i <- i+1
- if ((x >= X[i]-xinc) & (x < X[i]+xinc))
- break
- }
- ip <- i
- return(ip)
- }
-
Modified: pkg/R/is.3dpoints.R
===================================================================
--- pkg/R/is.3dpoints.R 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/is.3dpoints.R 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,9 +1,9 @@
-is.3dpoints <- function (p)
+is.3dpoints <- function (x)
{
is <- FALSE
- if (is.array(p))
- if (length(dim(p)) == 2)
- if (dim(p)[2] >= 3)
+ if (is.array(x))
+ if (length(dim(x)) == 2)
+ if (dim(x)[2] >= 3)
is <- TRUE
is
}
Deleted: pkg/R/make.grid.r
===================================================================
--- pkg/R/make.grid.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/make.grid.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,81 +0,0 @@
-make.grid <- function(nx,ny,poly)
- {
- #
- # Generate a rectangular grid of points in a polygon
- #
- # Requires Splancs package.
- #
- # Arguments:
- #
- # nx, ny: grid dimensions.
- #
- # poly: two columns matrix specifying polygonal region containing
- # all data locations. If poly is missing, the unit square is
- # considered.
- #
- # Value:
- #
- # x, y: numeric vectors giving the coordinates of the points
- # of the rectangular grid,
- #
- # X, Y: matrix containing x and y coordinates respectively,
- #
- # pts: index of the grid points belonging to the polygon,
- #
- # xinc, yinc: mesh size of the grid,
- #
- # mask: nx*ny matrix of logicals. TRUE if the grid point belongs
- # to the polygon.
- #
- ##
- ## E. GABRIEL, 28/12/2005
- ##
-
-# library(splancs)
-
- if (missing(poly)) poly <- matrix(c(0,0,1,0,1,1,0,1),4,2,T)
-
- if ((nx < 2) || (ny < 2)) stop("the grid must be at least of size 2x2")
-
- xrang <- range(poly[, 1], na.rm = TRUE)
- yrang <- range(poly[, 2], na.rm = TRUE)
- xmin <- xrang[1]
- xmax <- xrang[2]
- ymin <- yrang[1]
- ymax <- yrang[2]
-
- xinc <- (xmax-xmin)/nx
- yinc <- (ymax-ymin)/ny
-
- xc <- xmin-xinc/2
- yc <- ymin-yinc/2
- xgrid <- rep(0,nx)
- ygrid <- rep(0,ny)
- xgrid[1] <- xc + xinc
- ygrid[1] <- yc + yinc
-
- for (i in 2:nx)
- {
- xgrid[i] <- xgrid[i-1]+xinc
- }
- for (i in 2:ny)
- {
- ygrid[i] <- ygrid[i-1]+yinc
- }
-
- yy <- matrix(xgrid,nx,ny)
- xx <- t(yy)
- yy <- matrix(ygrid,nx,ny)
-
- X <- as.vector(xx)
- Y <- as.vector(yy)
-
- poly <- rbind(poly,poly[1,])
- pts <- inpip(pts=cbind(X,Y),poly)
-
- X[pts] <- TRUE
- X[X!=TRUE] <- FALSE
- mask <- matrix(X,ncol=ny,nrow=nx,byrow=T)
-
- invisible(return(list(x=xgrid,y=ygrid,X=xx,Y=yy,pts=pts,xinc=xinc,yinc=yinc,mask=matrix(as.logical(mask),nx,ny))))
- }
Deleted: pkg/R/pcp.larger.region.r
===================================================================
--- pkg/R/pcp.larger.region.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/pcp.larger.region.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,196 +0,0 @@
-pcp.larger.region <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infecD=TRUE, maxradlarger, ...)
-{
- if (missing(cluster)) cluster <- "uniform"
-
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
- if (missing(t.region)) t.region <- c(0,1)
-
- if (missing(maxrad)) maxrad <- c(0.05,0.05)
- if (length(maxrad)==1) maxrad=rep(maxrad,2)
- maxrads <- maxrad[1]
- maxradt <- maxrad[2]
-
- if (missing(maxradlarger)) maxradlarger <- maxrad
- maxradlargers <- maxradlarger[1]
- maxradlargert <- maxradlarger[2]
-
- if (length(cluster)==1)
- {
- s.distr <- cluster
- t.distr <- cluster
- }
- else
- {
- s.distr <- cluster[1]
- t.distr <- cluster[2]
- }
-
- t.region <- sort(t.region)
- s.area <- areapl(s.region)
- t.area <- t.region[2]-t.region[1]
-
- if (maxradlargers==0)
- s.larger <- s.region
- else
- {
- s.larger <- rbind(s.region+maxradlargers,s.region-maxradlargers,cbind(s.region[,1]+maxradlargers,s.region[,2]-maxradlargers),cbind(s.region[,1]-maxradlargers,s.region[,2]+maxradlargers))
- M <- chull(s.larger)
- s.larger <- cbind(s.larger[M,1],s.larger[M,2])
- }
- t.larger <- c(t.region[1]-maxradlargert,t.region[2]+maxradlargert)
-
- if (t.larger[1]<0) t.larger[1] <- 0
-
- if (is.null(nparents))
- {
- if (is.null(lambda)) stop("please specify either the number of parents or the intensity of parents process")
- if (is.function(lambda)) stop("please specify the number of parents")
- if (is.numeric(lambda))
- lambda <- lambda*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
- npar <- rpois(1,lambda)
- }
- else npar <- nparents
-
- if (is.null(lambda)) lambda <- (npar/(s.area*t.area))*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
-
- if (is.null(npoints))
- {
- if (is.null(mc)) stop("please specify either the number of points to simulate or the mean number of children per parents")
- else
- {
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- npts <- sum(nchild)
- }
- }
- else
- {
- npts <- npoints
- if (is.null(mc)) mc <- npts/npar
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- nc <- sum(nchild)
- while (nc < npts)
- {
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- nc <- sum(nchild)
- }
- }
-
- parpts <- rpp(lambda=lambda,s.region=s.larger,t.region=t.larger,npoints=npar,...)$xyt
-
- pattern.interm <- NULL
- ipt <- 1
- nchild2 <- rep(0,npar)
- for (ipar in 1:npar)
- {
- xpar <- parpts[ipar,1]
- ypar <- parpts[ipar,2]
- zpar <- parpts[ipar,3]
- nc <- nchild[ipar,1]
-
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(nc,min=-maxrads,max=maxrads)
- yp <- ypar+runif(nc,min=-maxrads,max=maxrads)
- }
- if (s.distr=="normal")
- {
- xp <- xpar+rnorm(nc,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(nc,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+rexp(nc,rate=1/maxrads)
- yp <- ypar+rexp(nc,rate=1/maxrads)
- }
- if (t.distr=="uniform")
- {
- if (infecD==TRUE)
- zp <- zpar+runif(nc,min=-maxradt,max=maxradt)
- else
- zp <- zpar+runif(nc,min=0,max=maxradt)
- }
- if (t.distr=="normal")
- {
- if (infecD==TRUE)
- zp <- zpar+abs(rnorm(nc,mean=0,sd=maxradt/2))
- else
- zp <- zpar+rnorm(nc,mean=0,sd=maxradt/2)
- }
- if (t.distr=="exponential")
- {
- if (infecD==TRUE)
- zp <- zpar+rexp(nc,rate=1/maxradt)
- else
- zp <- zpar+sample(c(-1,1),1)*rexp(nc,rate=1/maxradt)
- }
-
- mask <- ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
- nchild2[ipar] <- sum(mask)
- pattern.interm <- rbind(pattern.interm,cbind(xp[mask],yp[mask],zp[mask],rep(ipar,sum(mask))))
- }
-
- ipt <- dim(pattern.interm)[[1]]
- if (ipt > npts)
- {
- samp <- sample(1:ipt,npts)
- pattern.interm <- pattern.interm[samp,]
- }
- while(ipt < npts)
- {
- ipar <- sample(1:npar,1)
- nchild2[ipar] <- nchild2[ipar]+1
- xpar <- parpts[ipar,1]
- ypar <- parpts[ipar,2]
- zpar <- parpts[ipar,3]
-
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(1,min=-maxrads,max=maxrads)
- yp <- ypar+runif(1,min=-maxrads,max=maxrads)
- }
- if (s.distr=="normal")
- {
- xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+rexp(1,rate=1/maxrads)
- yp <- ypar+rexp(1,rate=1/maxrads)
- }
- if (t.distr=="uniform")
- {
- if (infecD==TRUE)
- zp <- zpar+runif(1,min=-maxradt,max=maxradt)
- else
- zp <- zpar+runif(1,min=0,max=maxradt)
- }
- if (t.distr=="normal")
- {
- if (infecD==TRUE)
- zp <- zpar+abs(rnorm(1,mean=0,sd=maxradt/2))
- else
- zp <- zpar+rnorm(1,mean=0,sd=maxradt/2)
- }
- if (t.distr=="exponential")
- {
- if (infecD==TRUE)
- zp <- zpar+rexp(1,rate=1/maxradt)
- else
- zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/maxradt)
- }
-
- if ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
- {
- pattern.interm <- rbind(pattern.interm,cbind(xp,yp,zp,ipar))
- ipt <- ipt+1
- }
- }
-
- pts <- pattern.interm[,1:3]
- ott<-order(pts[,3])
- pts<-pts[ott,]
- invisible(return(list(pts=pts,s.larger=s.larger,t.larger=t.larger,index.child=pattern.interm[,4],nchild=nchild2,parpts=parpts)))
-}
-
-
Modified: pkg/R/plot.stpp.R
===================================================================
--- pkg/R/plot.stpp.R 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/plot.stpp.R 2010-02-13 13:49:51 UTC (rev 31)
@@ -5,11 +5,11 @@
{
par(mfrow=c(1,2),pty="s")
if (is.null(s.region))
- plot(x[,1:2],pch=20,main="xy-locations")
+ plot(x[,1:2],main="xy-locations",...)
else
{
polymap(s.region,xlab="x",ylab="y")
- points(x[,1:2],pch=20)
+ points(x[,1:2],...)
title("xy-locations")
}
plot(x[,3],cumsum(x[,3]),type="l",xlab="t",ylab="",main="cumulative number",las=1,xlim=t.region)
Deleted: pkg/R/rhpp.r
===================================================================
--- pkg/R/rhpp.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/rhpp.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,59 +0,0 @@
-rhpp <- function(lambda, s.region, t.region, npoints=NULL, replace=TRUE, discrete.time=FALSE)
- {
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
- if (missing(t.region)) t.region <- c(0,1)
-
- s.area <- areapl(s.region)
- t.region <- sort(t.region)
- t.area <- t.region[2]-t.region[1]
-
- if (missing(lambda) & !(is.null(npoints)))
- {
- if (t.area==0) lambda <- npoints/(s.area)
- else lambda <- npoints/(s.area * t.area)
- }
-
- pattern <- list()
- index.t <- list()
- if (is.numeric(lambda))
- {
- if (is.null(npoints)==TRUE)
- {
- if (t.area==0)
- { npoints <- round(rpois(n=1,lambda=lambda * s.area),0) }
- else
- { npoints <- round(rpois(n=1,lambda=lambda * s.area * t.area),0) }
- }
- xy <- matrix(csr(poly=s.region,npoints=npoints),ncol=2)
- x <- xy[,1]
- y <- xy[,2]
- npoints <- length(x)
-
- if (discrete.time==TRUE)
- {
- vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
- if ((length(vect)<npoints) & (replace==F))
- stop("when replace=FALSE and discrete.time=TRUE, the length of seq(t.region[1],t.region[2],by=1) must be greater than the number of points")
- names(vect) <- 1:length(vect)
- M <- sample(vect,npoints,replace=replace)
- times <- M
- names(times) <- NULL
- samp <- as.numeric(names(M))
- }
- else
- {
- times <- runif(npoints,min=t.region[1],max=t.region[2])
- samp <- sample(1:npoints,npoints,replace=replace)
- times <- times[samp]
- }
- times <- sort(times)
- index.times <- sort(samp)
- pattern.interm <- list(x=x,y=y,t=times,index.t=index.times)
- pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)
- index.t <- pattern.interm$index.t
- }
- else
- stop("lambda must be numeric")
-
- invisible(return(list(pts=pattern,index.t=index.t)))
- }
Modified: pkg/R/rinter.r
===================================================================
--- pkg/R/rinter.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/rinter.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,3 +1,358 @@
+spatial.inhibition <- function(npoints,h,theta,delta,p,recent="all",s.region,inhibition=TRUE,...)
+ {
+ #
+ # Simulate an inhibition or a contagious spatial point process in a region D
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # s.region: two columns matrix specifying polygonal region containing
+ # all data locations. If poly is missing, the unit square is
+ # considered.
+ #
+ # h: a function of the distance between points.
+ # If inhibition=TRUE, h is monotone, increasing, and must tend
+ # to 1 when the distance tends to infinity. 0<= h(d,theta,delta) <= 1.
+ # Else, h is monotone, decreasing, and must tend
+ # to 1 when the distance tends to 0. 0<= h(d,theta,delta) <= 1.
+ #
+ # p: a function among "min", "max", "prod".
+ #
+ # npoints: number of points to simulate.
+ #
+ #
+ # Value:
+ # pattern: list containing the points (x,y) of the simulated point
+ # process.
+ #
+ ##
+ ## E. GABRIEL, 26/03/2007
+ ##
+ ## last modification: 21/10/2008
+ ##
+ ##
+
+ if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
+
+ if (!(is.function(h)))
+ {
+ models <- c("step","gaussian")
+ if (sum(h==models)==0)
+ {
+ message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
+ stop(message)
+ }
+ if (h=="step")
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- rep(1,length(d))
+ if (inhibition==TRUE) res[d<=delta] <- theta
+ else res[d>=delta] <- theta
+ return(res)
+ }
+ }
+ if (h=="gaussian")
+ {
+ hk <- function(d,theta,delta)
+ {
+ if (inhibition==TRUE)
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,0)
+ if (d[i]>(delta+theta/2)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+theta/2)) res=c(res,exp(-((d[i]-delta-theta/2)^2)/(2*(theta/8)^2)))
+ }
+ }
+ else
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<delta) res=c(res,1)
+ else res=c(res,exp(-((d[i]-delta)^2)/(2*(theta/8)^2)))
+ }
+ }
+ return(res)
+ }
+ }
+ }
+ else
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- h(d,theta,delta)
+ return(res)
+ }
+ }
+
+ pk <- function(d,h,recent,theta,delta)
+ {
+ if (recent=="all")
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(d))
+ {
+ if (p=="min") res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ }
+ else
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ }
+ else
+ stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+
+ xy <- csr(n=1,poly=s.region)
+ npts <- 1
+ pattern.interm <- cbind(x=xy[1],y=xy[2])
+ if (inhibition==TRUE)
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+ xy <- csr(n=1,poly=s.region)
+ if (all((sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2)) > delta))
+ umax <- 1
+ else
+ umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),hk,recent,theta,delta)
+
+ if (prob<umax)
+ {
+ pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
+ npts <- npts+1
+ }
+ }
+ }
+ else
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+ continue <- FALSE
+ while(continue==FALSE)
+ {
+ xy <- csr(n=1,poly=s.region)
+ if (all(sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2) < delta))
+ umax <- 1
+ else
+ umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),hk,recent,theta,delta)
+
+ if (prob < umax)
+ {
+ pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
+ npts <- npts+1
+ continue <- TRUE
+ }
+ }
+ }
+ }
+ invisible(return(list(pts=pattern.interm,s.region=s.region)))
+}
+
+
+temporal.inhibition <- function(npoints,h,theta,delta,p,recent="all",t.region,discrete.time=FALSE,replace=FALSE,inhibition=TRUE)
+ {
+ #
+ # Simulate an inhibition or a contagious temporal point process in T
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # t.region: vector containing the minimum and maximum values of
+ # the time interval.
+ #
+ # h: a function of the distance between times and theta.
+ # If inhibition=TRUE, h is monotone, increasing, and must tend
+ # to 1 when the distance tends to infinity. 0<= h(d,theta,delta) <= 1.
+ # Else, h is monotone, decreasing, and must tend
+ # to 1 when the distance tends to 0. 0<= h(d,theta,delta) <= 1.
+ #
+ # p: a function among "min", "max", "prod".
+ #
+ # replace: logical allowing times repetition.
+ #
+ # npoints: number of points to simulate.
+ #
+ #
+ # Value:
+ # pattern: list containing times t of the simulated process.
+ #
+ ##
+ ## E. GABRIEL, 26/03/2007
+ ##
+ ## last modification: 31/10/2008
+ ##
+ ##
+
+ if (missing(t.region)) t.region <- c(0,1)
+
+ if (!(is.function(h)))
+ {
+ models <- c("step","gaussian")
+ if (sum(h==models)==0)
+ {
+ message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
+ stop(message)
+ }
+ if (h=="step")
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- rep(1,length(d))
+ if (inhibition==TRUE) res[d<=delta] <- theta
+ else res[d>=delta] <- theta
+ return(res)
+ }
+ }
+ if (h=="gaussian")
+ {
+ hk <- function(d,theta,delta)
+ {
+ if (inhibition==TRUE)
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,0)
+ if (d[i]>(delta+theta/2)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+theta/2)) res=c(res,exp(-((d[i]-delta-theta/2)^2)/(2*(theta/8)^2)))
+ }
+ }
+ else
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<delta) res=c(res,1)
+ else res=c(res,exp(-((d[i]-delta)^2)/(2*(theta/8)^2)))
+ }
+ }
+ return(res)
+ }
+ }
+ }
+ else
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- h(d,theta,delta)
+ return(res)
+ }
+ }
+
+ pk <- function(d,h,recent,theta,delta)
+ {
+ if (recent=="all")
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(d))
+ {
+ if (p=="min") res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ }
+ else
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ }
+ else stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[1]+delta)
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[1]+delta),1)
+ times <- ti
+ npts <- 1
+ if (inhibition==TRUE)
+ {
+ while(npts < npoints)
+ {
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[2])
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
+
+ prob <- runif(1)
+
+ if (all(abs(ti - times) > delta))
+ umax <- 1
+ else
+ umax <- pk(d=abs(ti - times),hk,recent,theta,delta)
+ if (prob<umax)
+ {
+ times <- c(times,ti)
+ npts <- npts+1
+ }
+ }
+ }
+ else
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+
+ continue <- FALSE
+ while(continue==FALSE)
+ {
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[2])
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
+
+ if (abs(ti - times[npts]) < delta)
+ umax <- 1
+ else
+ umax <- pk(d=abs(ti - times),hk,recent,theta,delta)
+ if (prob < umax)
+ {
+ times <- c(times,ti)
+ npts <- npts+1
+ continue <- TRUE
+ }
+ }
+ }
+ }
+
+ samp <- sample(1:npoints,npoints,replace=replace)
+ times <- sort(times[samp])
+
+ invisible(return(list(times=times,t.region=t.region)))
+}
+
+
+
rinter <- function(npoints,s.region,t.region,hs="step",gs="min",thetas=0,deltas,ht="step",gt="min",thetat=1,deltat,recent="all",nsim=1,discrete.time=FALSE,replace=FALSE,inhibition=TRUE,...)
{
#
Deleted: pkg/R/ripp.r
===================================================================
--- pkg/R/ripp.r 2010-02-12 16:33:59 UTC (rev 30)
+++ pkg/R/ripp.r 2010-02-13 13:49:51 UTC (rev 31)
@@ -1,286 +0,0 @@
-ripp <- function(lambda, s.region, t.region, npoints=NULL, replace=T, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=NULL, mut=NULL, ...)
-{
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
- if (missing(t.region)) t.region <- c(0,1)
-
- s.area <- areapl(s.region)
- t.region <- sort(t.region)
- t.area <- t.region[2]-t.region[1]
-
- if (missing(lambda) & !(is.null(npoints)))
- {
- if (t.area==0) lambda <- npoints/(s.area)
- else lambda <- npoints/(s.area * t.area)
- }
-
- lambdamax <- lmax
- pattern <- list()
- index.t <- list()
-
- if (is.function(lambda))
- {
- s.grid <- make.grid(nx,ny,s.region)
- s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
- if (discrete.time==T)
- {
- vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
- if (nt>length(vect))
- {
- nt <- length(vect)
- warning("nt used is less than the one given in argument")
- t.grid <- list(times=vect,tinc=1)
- }
- else
- {
- vect <- round(seq(floor(t.region[1]),ceiling(t.region[2]),length=nt))
- t.grid <- list(times=vect,tinc=round(t.area/(nt-1)))
- }
- }
- else
- t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
-
- if (is.null(Lambda) & is.null(mut))
- {
- Lambda <- array(NaN,dim=c(nx,ny,nt))
- mut <- rep(0,nt)
- # maxlambda <- rep(0,nt)
- for(it in 1:nt)
- {
- L <- lambda(as.vector(s.grid$X),as.vector(s.grid$Y),t.grid$times[it],...)
- M <- matrix(L,ncol=ny,nrow=nx,byrow=T)
- M[!(s.grid$mask)] <- NaN
- Lambda[,,it] <- M
- # maxlambda[it] <- max(Lambda[,,it],na.rm=T)
- mut[it] <- sum(Lambda[,,it],na.rm=T)
- }
- }
-
- if (is.null(npoints)==T)
- {
- if (t.area==0)
- { en <- sum(Lambda,na.rm=T)*s.grid$xinc*s.grid$yinc }
- else
- {
- en <- sum(Lambda,na.rm=T)*s.grid$xinc*s.grid$yinc*t.grid$tinc
- npoints <- round(rpois(n=1,lambda=en),0)
- }
- }
-
- if (is.null(lambdamax))
- lambdamax <- max(Lambda,na.rm=T)
- npts <- round(lambdamax/(s.area*t.area),0)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/stpp -r 31
More information about the Stpp-commits
mailing list