[Stpp-commits] r60 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 21 11:25:25 CEST 2012
Author: gabriele
Date: 2012-08-21 11:25:25 +0200 (Tue, 21 Aug 2012)
New Revision: 60
Modified:
pkg/DESCRIPTION
pkg/R/rinfec.r
pkg/R/rlgcp.r
pkg/R/rpp.r
pkg/man/rinfec.Rd
pkg/man/stpp-package.Rd
Log:
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-08-20 08:02:56 UTC (rev 59)
+++ pkg/DESCRIPTION 2012-08-21 09:25:25 UTC (rev 60)
@@ -1,13 +1,13 @@
-Package: stpp
-Type: Package
-Title: Space-Time Point Pattern simulation, visualisation and analysis
-Version: 1.0
-Date: 2012-08-08
-Author: Edith Gabriel, Peter J Diggle, stan function by Barry Rowlingson
-Maintainer: Edith Gabriel <edith.gabriel at univ-avignon.fr>
-Depends: R (>= 2.10), splancs, KernSmooth
-Suggests: rpanel, rgl
-Description: A package for analysing, simulating and displaying space-time point patterns
-License: GPL-3
-LazyLoad: yes
-
+Package: stpp
+Type: Package
+Title: Space-Time Point Pattern simulation, visualisation and analysis
+Version: 1.0-1
+Date: 2012-08-08
+Author: Edith Gabriel, Peter J Diggle, stan function by Barry Rowlingson
+Maintainer: Edith Gabriel <edith.gabriel at univ-avignon.fr>
+Depends: R (>= 2.10), splancs, KernSmooth
+Suggests: rpanel, rgl
+Description: A package for analysing, simulating and displaying space-time point patterns
+License: GPL-3
+LazyLoad: yes
+
Modified: pkg/R/rinfec.r
===================================================================
--- pkg/R/rinfec.r 2012-08-20 08:02:56 UTC (rev 59)
+++ pkg/R/rinfec.r 2012-08-21 09:25:25 UTC (rev 60)
@@ -226,13 +226,17 @@
up <- runif(1)
nix <- findInterval(vec=s.grid$x,x=xp)
niy <- findInterval(vec=s.grid$y,x=yp)
- Up <- Lambda[nix,niy]/lmax
- retain <- up <= Up
- if ((retain==FALSE) || is.na(retain))
- retain.eq.F <- FALSE
- else
- retain.eq.F <- TRUE
-
+ if (nix!=0 & niy!=0)
+ {
+ Up <- Lambda[nix,niy]/lmax
+ retain <- up <= Up
+ if ((retain==FALSE) || is.na(retain))
+ retain.eq.F <- FALSE
+ else
+ retain.eq.F <- TRUE
+ }
+ else
+ retain.eq.F <- FALSE
}
}
M <- inout(pts=rbind(c(xp,yp),c(xp,yp)),poly=s.region)
Modified: pkg/R/rlgcp.r
===================================================================
--- pkg/R/rlgcp.r 2012-08-20 08:02:56 UTC (rev 59)
+++ pkg/R/rlgcp.r 2012-08-21 09:25:25 UTC (rev 60)
@@ -354,7 +354,10 @@
nix <- findInterval(vec=s.grid$x,x=x[ix])
niy <- findInterval(vec=s.grid$y,x=y[ix])
nit <- findInterval(vec=t.grid$times,x=times[ix])
- prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
+ if (nix==0 | niy==0 | nit==0)
+ prob=c(prob,NA)
+ else
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
@@ -407,7 +410,10 @@
nix <- findInterval(vec=s.grid$x,x=wx[ix])
niy <- findInterval(vec=s.grid$y,x=wy[ix])
nit <- findInterval(vec=t.grid$times,x=wtimes[ix])
- prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
+ if (nix==0 | niy==0 | nit==0)
+ prob=c(prob,NA)
+ else
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
if (length(M)!=0)
Modified: pkg/R/rpp.r
===================================================================
--- pkg/R/rpp.r 2012-08-20 08:02:56 UTC (rev 59)
+++ pkg/R/rpp.r 2012-08-21 09:25:25 UTC (rev 60)
@@ -1,534 +1,540 @@
-.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)
-
- 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==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)
-
- 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==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)
- }
- 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
- }
-
- 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)
-# npts <- round(lambdamax/(s.area*t.area),0)
- npts <- npoints
- 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 (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]
-
- retain.eq.F <- FALSE
- while(retain.eq.F==FALSE)
- {
- xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
- x <- xy[,1]
- y <- xy[,2]
-
- prob <- NULL
- for(nx in 1:length(x))
- {
- nix <- findInterval(vec=s.grid$x,x=x[nx])
- niy <- findInterval(vec=s.grid$y,x=y[nx])
- nit <- findInterval(vec=t.grid$times,x=times[nx])
- prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
- }
-
- M <- which(is.na(prob))
- if (length(M)!=0)
- {
- x <- x[-M]
- y <- y[-M]
- times <- times[-M]
- prob <- prob[-M]
- npts <- length(x)
- }
-
- u <- runif(npts)
- retain <- u <= prob
- if (sum(retain==F)==length(retain)) retain.eq.F <- FALSE
- else retain.eq.F <- TRUE
- }
-# if (sum(retain==FALSE)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
-
-
- 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==FALSE)==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)
- if (neffec > npoints)
- {
- retain <- 1:npoints
- x <- x[retain]
- y <- y[retain]
- samp <- samp[retain]
- samp.remain <- (1:nt)[-samp]
- times <- times[retain]
- }
-
- 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 <- NULL
- for(nx in 1:length(wx))
- {
- nix <- findInterval(vec=s.grid$x,x=wx[nx])
- niy <- findInterval(vec=s.grid$y,x=wy[nx])
- nit <- findInterval(vec=t.grid$times,x=wtimes[nx])
- prob <- c(prob,(Lambda[nix,niy,nit])/lambdamax)
- }
- M <- which(is.na(prob))
- if (length(M)!=0)
- {
- wx <- wx[-M]
- wy <- wy[-M]
- wtimes <- wtimes[-M]
- prob <- prob[-M]
- }
- if (neffec > 0)
- {
- u <- runif(length(prob))
- 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)
- }
- }
- 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
- }
-
- invisible(return(list(pts=pattern,index.t=index.t)))
-}
-
-
-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)
-
- 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()
- 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)
- index.t <- ipp$index.t
- }
- else
- {
- pattern[[ni]] <- as.3dpoints(ipp$pts)
- index.t[[ni]] <- ipp$index.t
- }
- 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)
- index.t <- ipp$index.t
- }
- else
- {
- pattern[[ni]] <- as.3dpoints(ipp$pts)
- index.t[[ni]] <- ipp$index.t
- }
- 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,index.t=index.t,s.region=s.region,t.region=t.region,lambda=lambda,Lambda=Lambda)))
-}
-
-
-
+.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)
+
+ 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==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)
+
+ 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==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)
+ }
+ 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
+ }
+
+ 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)
+# npts <- round(lambdamax/(s.area*t.area),0)
+ npts <- npoints
+ 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 (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]
+
+ retain.eq.F <- FALSE
+ while(retain.eq.F==FALSE)
+ {
+ xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
+ x <- xy[,1]
+ y <- xy[,2]
+
+ prob <- NULL
+ for(nx in 1:length(x))
+ {
+ nix <- findInterval(vec=s.grid$x,x=x[nx])
+ niy <- findInterval(vec=s.grid$y,x=y[nx])
+ nit <- findInterval(vec=t.grid$times,x=times[nx])
+ if (nix==0 | niy==0 | nit==0)
+ prob=c(prob,NA)
+ else
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
+ }
+
+ M <- which(is.na(prob))
+ if (length(M)!=0)
+ {
+ x <- x[-M]
+ y <- y[-M]
+ times <- times[-M]
+ prob <- prob[-M]
+ npts <- length(x)
+ }
+
+ u <- runif(npts)
+ retain <- u <= prob
+ if (sum(retain==F)==length(retain)) retain.eq.F <- FALSE
+ else retain.eq.F <- TRUE
+ }
+# if (sum(retain==FALSE)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
+
+
+ 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==FALSE)==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)
+ if (neffec > npoints)
+ {
+ retain <- 1:npoints
+ x <- x[retain]
+ y <- y[retain]
+ samp <- samp[retain]
+ samp.remain <- (1:nt)[-samp]
+ times <- times[retain]
+ }
+
+ 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]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/stpp -r 60
More information about the Stpp-commits
mailing list