Nothing
.make.grid <- function(nx,ny,poly)
{
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=TRUE)
invisible(return(list(x=xgrid,y=ygrid,X=xx,Y=yy,pts=pts,xinc=xinc,yinc=yinc,mask=matrix(as.logical(mask),nx,ny))))
}
.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)
xp <- s.region[,1]
yp <- s.region[,2]
nedges <- length(xp)
yp <- yp - min(yp)
nxt <- c(2:nedges, 1)
dx <- xp[nxt] - xp
ym <- (yp + yp[nxt])/2
Areaxy <- -sum(dx * ym)
if (Areaxy > 0){
bdry <- owin(poly = list(x = s.region[,1], y = s.region[,2]))
}else{
bdry <- owin(poly = list(x = s.region[,1][length(s.region[,1]):1], y = s.region[,2][length(s.region[,1]):1]))
}
s.area <- area(bdry)
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==FALSE))
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)))
}
.ripp <- function(lambda, s.region, t.region, npoints=NULL, replace=TRUE, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=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)
xp <- s.region[,1]
yp <- s.region[,2]
nedges <- length(xp)
yp <- yp - min(yp)
nxt <- c(2:nedges, 1)
dx <- xp[nxt] - xp
ym <- (yp + yp[nxt])/2
Areaxy <- -sum(dx * ym)
if (Areaxy > 0){
bdry <- owin(poly = list(x = s.region[,1], y = s.region[,2]))
}else{
bdry <- owin(poly = list(x = s.region[,1][length(s.region[,1]):1], y = s.region[,2][length(s.region[,1]):1]))
}
s.area <- area(bdry)
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==TRUE)
{
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))
{
Lambda <- array(NaN,dim=c(nx,ny,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=TRUE)
M[!(s.grid$mask)] <- NaN
Lambda[,,it] <- M
}
}
if (is.null(npoints)==TRUE)
{
if (t.area==0)
{ en <- sum(Lambda,na.rm=TRUE)*s.grid$xinc*s.grid$yinc }
else
{
en <- sum(Lambda,na.rm=TRUE)*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=TRUE)
npts <- round(lambdamax/(s.area*t.area),0)
if (npts==0) stop("there is no data to thin")
xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
x <- xy[,1]
y <- xy[,2]
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)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
times.init <- sample(vect,nt,replace=replace)
}
else
times.init <- runif(nt,min=t.region[1],max=t.region[2])
samp <- sample(1:nt,npts,replace=replace)
times <- times.init[samp]
prob <- lambda(x,y,times,...)/lambdamax
u <- runif(npts)
retain <- u <= prob
if (sum(retain==FALSE)==length(retain))
{
lambdas <- matrix(0,nrow=nx,ncol=ny)
for(ix in 1:nx){for(iy in 1:ny){
lambdas[ix,iy] <- median(Lambda[ix,iy,],na.rm=TRUE)}}
lambdamax <- max(lambdas,na.rm=TRUE)
prob <- lambda(x,y,times,...)/lambdamax
retain <- u <= prob
if (sum(retain==F)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
}
x <- x[retain]
y <- y[retain]
samp <- samp[retain]
samp.remain <- (1:nt)[-samp]
times <- times[retain]
neffec <- length(x)
while(neffec < npoints)
{
xy <- as.matrix(csr(poly=s.region,npoints=npoints-neffec))
if(dim(xy)[2]==1){wx <- xy[1]; wy <- xy[2]}
else{wx <- xy[,1]; wy <- xy[,2]}
if(replace==FALSE)
{ wsamp <- sample(samp.remain,npoints-neffec,replace=replace) }
else{ wsamp <- sample(1:nt,npoints-neffec,replace=replace) }
wtimes <- times.init[wsamp]
# lambdamax <- maxlambda[wsamp]
prob <- lambda(wx,wy,wtimes,...)/lambdamax
u <- runif(npoints-neffec)
retain <- u <= prob
x <- c(x,wx[retain])
y <- c(y,wy[retain])
times <- c(times,wtimes[retain])
samp <- c(samp,wsamp[retain])
samp.remain <- (1:nt)[-samp]
neffec <- length(x)
}
index.times <- sort(times,index.return=TRUE)$ix
pattern.interm <- list(x=x[index.times],y=y[index.times],t=times[index.times])
pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)[index.times,]
}
if (is.character(lambda))
{
if (is.null(Lambda))
stop("Lambda must be specified")
nx <- dim(Lambda)[1]
ny <- dim(Lambda)[2]
nt <- dim(Lambda)[3]
s.grid <- .make.grid(nx,ny,s.region)
s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
if (discrete.time==TRUE)
{
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(npoints))
{
en <- sum(Lambda,na.rm=TRUE)*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=TRUE)
if ((replace==FALSE) & (nt < npoints)) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==TRUE)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
times.init <- sample(vect,nt,replace=replace)
t.grid$times = times.init
}
### else
### times.init <- runif(nt,min=t.region[1],max=t.region[2])
XX=rep(s.grid$X,nt)
YY=rep(s.grid$Y,nt)
TT=rep(t.grid$times,each=length(s.grid$X))
df=NULL
for(nl in 1:nt)
{
LL = Lambda
LL[is.na(LL)] = -999
lambdal=as.im(list(x=s.grid$x,y=s.grid$y,z=LL[,,nl]))
df <- rbind(df,as.data.frame(lambdal))
}
df$value[df$value==-999]=0
samp=sample.int(length(XX),npoints,replace=TRUE,prob=df$value)
xx <- XX[samp] + runif(npoints, -s.grid$xinc/2, s.grid$xinc/2)
yy <- YY[samp] + runif(npoints, -s.grid$yinc/2, s.grid$yinc/2)
if (discrete.time==TRUE)
tt <- floor(TT[samp] + runif(npoints, -t.grid$tinc/2, t.grid$tinc/2))
else
tt <- TT[samp] + runif(npoints, -t.grid$tinc/2, t.grid$tinc/2)
xyt.init=cbind(x=xx,y=yy,t=tt)
retain.eq.F <- FALSE
while(retain.eq.F==FALSE)
{
pts <- inpip(pts=cbind(xyt.init[,1],xyt.init[,2]),poly=s.region)
xyt.init <- xyt.init[pts,]
ptt <- xyt.init[,3]>=t.region[1] & xyt.init[,3]<=t.region[2]
xyt.init <- xyt.init[ptt,]
npts = dim(xyt.init)[1]
if (npts == npoints) retain.eq.F <- TRUE
else
{
samp=sample.int(length(XX),npoints-npts,replace=TRUE,prob=df$value)
xx <- XX[samp] + runif(npoints-npts, -s.grid$xinc/2, s.grid$xinc/2)
yy <- YY[samp] + runif(npoints-npts, -s.grid$yinc/2, s.grid$yinc/2)
if (discrete.time==TRUE)
tt <- floor(TT[samp] + runif(npoints-npts, -t.grid$tinc/2, t.grid$tinc/2))
else
tt <- TT[samp] + runif(npoints-npts, -t.grid$tinc/2, t.grid$tinc/2)
xyt.init1=cbind(x=xx,y=yy,t=tt)
xyt.init=rbind(xyt.init,xyt.init1)
}}
x<-xyt.init[,1]
y<-xyt.init[,2]
times<-xyt.init[,3]
index.times <- sort(times,index.return=TRUE)$ix
pattern.interm <- list(x=x[index.times],y=y[index.times],t=times[index.times])
pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)[index.times,]
}
invisible(return(list(pts=pattern)))
}
rpp <- function(lambda, s.region, t.region, npoints=NULL, nsim=1, replace=TRUE, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=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)
xp <- s.region[,1]
yp <- s.region[,2]
nedges <- length(xp)
yp <- yp - min(yp)
nxt <- c(2:nedges, 1)
dx <- xp[nxt] - xp
ym <- (yp + yp[nxt])/2
Areaxy <- -sum(dx * ym)
if (Areaxy > 0){
bdry <- owin(poly = list(x = s.region[,1], y = s.region[,2]))
}else{
bdry <- owin(poly = list(x = s.region[,1][length(s.region[,1]):1], y = s.region[,2][length(s.region[,1]):1]))
}
s.area <- area(bdry)
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()
ni <- 1
#
# Homogeneous Poisson Process
#
if (is.numeric(lambda) & length(lambda)==1)
{
while(ni<=nsim)
{
hpp <- .rhpp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time)
if (nsim==1)
{
pattern <- as.3dpoints(hpp$pts)
index.t <- hpp$index.t
}
else
{
pattern[[ni]] <- as.3dpoints(hpp$pts)
index.t[[ni]] <- hpp$index.t
}
ni <- ni+1
}
Lambda <- NULL
}
#
# Inhomogeneous Poisson Process
#
else 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==TRUE)
{
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)))
Lambda <- array(NaN,dim=c(nx,ny,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=TRUE)
M[!(s.grid$mask)] <- NaN
Lambda[,,it] <- M
}
while(ni<=nsim)
{
ipp <- .ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, ...)
if (nsim==1)
{
pattern <- as.3dpoints(ipp$pts)
}
else
{
pattern[[ni]] <- as.3dpoints(ipp$pts)
}
ni <- ni+1
}
}
else if (is.array(lambda))
{
if (length(dim(lambda))!=3) stop ("lambda must be a 3D-array")
Lambda = lambda
lambda = "a"
while(ni<=nsim)
{
ipp <- .ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, ...)
if (nsim==1)
{
pattern <- as.3dpoints(ipp$pts)
}
else
{
pattern[[ni]] <- as.3dpoints(ipp$pts)
}
ni <- ni+1
}
}
else stop("lambda must be either a single positive value or a function or a 3D-array")
invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region,lambda=lambda,Lambda=Lambda)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.