# nndist.R # Check that nndist and nnwhich give # results consistent with direct calculation from pairdist # Similarly for nncross and distfun require(spatstat) local({ eps <- sqrt(.Machine$double.eps) f <- function(mat,k) { apply(mat, 1, function(z,n) { sort(z)[n] }, n=k+1) } g <- function(mat,k) { apply(mat, 1, function(z,n) { order(z)[n] }, n=k+1) } # Two dimensions X <- runifpoint(42) nn <- nndist(X) nnP <- f(pairdist(X), 1) if(any(abs(nn - nnP) > eps)) stop("nndist.ppp does not agree with pairdist") nn5 <- nndist(X, k=5) nn5P <- f(pairdist(X), 5) if(any(abs(nn5 - nn5P) > eps)) stop("nndist.ppp(k=5) does not agree with pairdist") nw <- nnwhich(X) nwP <- g(pairdist(X), 1) if(any(nw != nwP)) stop("nnwhich.ppp does not agree with pairdist") nw5 <- nnwhich(X, k=5) nw5P <- g(pairdist(X), 5) if(any(nw5 != nw5P)) stop("nnwhich.ppp(k=5) does not agree with pairdist") X <- runifpoint3(42) nn <- nndist(X) nnP <- f(pairdist(X), 1) if(any(abs(nn - nnP) > eps)) stop("nndist.pp3 does not agree with pairdist") nn5 <- nndist(X, k=5) nn5P <- f(pairdist(X), 5) if(any(abs(nn5 - nn5P) > eps)) stop("nndist.pp3(k=5) does not agree with pairdist") nw <- nnwhich(X) nwP <- g(pairdist(X), 1) if(any(nw != nwP)) stop("nnwhich.pp3 does not agree with pairdist") nw5 <- nnwhich(X, k=5) nw5P <- g(pairdist(X), 5) if(any(nw5 != nw5P)) stop("nnwhich.pp3(k=5) does not agree with pairdist") # m dimensions X <- runifpointx(42, boxx(c(0,1),c(0,1),c(0,1),c(0,1))) nn <- nndist(X) nnP <- f(pairdist(X), 1) if(any(abs(nn - nnP) > eps)) stop("nndist.ppx does not agree with pairdist") nn5 <- nndist(X, k=5) nn5P <- f(pairdist(X), 5) if(any(abs(nn5 - nn5P) > eps)) stop("nndist.ppx(k=5) does not agree with pairdist") nw <- nnwhich(X) nwP <- g(pairdist(X), 1) if(any(nw != nwP)) stop("nnwhich.ppx does not agree with pairdist") nw5 <- nnwhich(X, k=5) nw5P <- g(pairdist(X), 5) if(any(nw5 != nw5P)) stop("nnwhich.ppx(k=5) does not agree with pairdist") #### nncross in two dimensions X <- runifpoint(42) Y <- runifpoint(42, win=owin(c(1,2),c(1,2))) # default nncross nc <- nncross(X,Y) ncd <- nc$dist ncw <- nc$which cd <- crossdist(X,Y) cdd <- apply(cd, 1, min) cdw <- apply(cd, 1, which.min) if(any(abs(ncd - cdd) > eps)) stop("nncross()$dist does not agree with apply(crossdist(), 1, min)") if(any(ncw != cdw)) stop("nncross()$which does not agree with apply(crossdist(), 1, which.min)") # sort on x nc <- nncross(X,Y, sortby="x") ncd <- nc$dist ncw <- nc$which if(any(abs(ncd - cdd) > eps)) stop("nncross(sortby=x)$dist does not agree with apply(crossdist(), 1, min)") if(any(ncw != cdw)) stop("nncross(sortby=x)$which does not agree with apply(crossdist(), 1, which.min)") # pre-sorted on x Y <- Y[order(Y$x)] nc <- nncross(X,Y, is.sorted.Y=TRUE, sortby="x") ncd <- nc$dist ncw <- nc$which cd <- crossdist(X,Y) cdd <- apply(cd, 1, min) cdw <- apply(cd, 1, which.min) if(any(abs(ncd - cdd) > eps)) stop("For sorted data, nncross()$dist does not agree with apply(crossdist(), 1, min)") if(any(ncw != cdw)) stop("For sorted data, nncross()$which does not agree with apply(crossdist(), 1, which.min)") }) require(spatstat) local({ Y <- split(urkiola) B <- Y$birch O <- Y$oak B.lam <- predict (ppm(B, ~polynom(x,y,2)), type="trend") O.lam <- predict (ppm(O, ~polynom(x,y,2)), type="trend") Kinhom(B, lambda=B.lam, correction="iso") Kinhom(B, lambda=B.lam, correction="border") Kcross.inhom(urkiola, i="birch", j="oak", B.lam, O.lam) Kcross.inhom(urkiola, i="birch", j="oak", B.lam, O.lam, correction = "iso") Kcross.inhom(urkiola, i="birch", j="oak", B.lam, O.lam, correction = "border") }) require(spatstat) local({ # critical R values that provoke GCC bug #323 a <- marktable(lansing, R=0.25) a <- marktable(lansing, R=0.21) a <- marktable(lansing, R=0.20) a <- marktable(lansing, R=0.10) }) # tests/NAinCov.R # Testing the response to the presence of NA's in covariates require(spatstat) local({ X <- runifpoint(42) Y <- as.im(function(x,y) { x+y }, owin()) Y[owin(c(0.2,0.4),c(0.2,0.4))] <- NA # fit model: should produce a warning but no failure misfit <- ppm(X, ~Y, covariates=list(Y=Y)) # prediction Z <- predict(misfit, type="trend") Z <- predict(misfit, type="se") # covariance matrix: all should be silent v <- vcov(misfit) ss <- vcov(misfit, what="internals") NULL }) # tests for agreement between C and interpreted code # for interpoint distances require(spatstat) local({ eps <- .Machine$double.eps * 4 # pairdist.ppp X <- rpoispp(42) dC <- pairdist(X, method="C") dR <- pairdist(X, method="interpreted") if(any(abs(dC - dR) > eps)) stop("Algorithms for pairdist() do not agree") dC <- pairdist(X, periodic=TRUE, method="C") dR <- pairdist(X, periodic=TRUE, method="interpreted") if(any(abs(dC - dR) > eps)) stop("Algorithms for pairdist(periodic=TRUE) do not agree") # crossdist.ppp Y <- rpoispp(42) dC <- crossdist(X, Y, method="C") dR <- crossdist(X, Y, method="interpreted") if(any(abs(dC - dR) > eps)) stop("Algorithms for crossdist() do not agree") dC <- crossdist(X, Y, periodic=TRUE, method="C") dR <- crossdist(X, Y, periodic=TRUE, method="interpreted") if(any(abs(dC - dR) > eps)) stop("Algorithms for crossdist(periodic=TRUE) do not agree") # nndist.ppp nnC <- nndist(X, method="C") nnI <- nndist(X, method="interpreted") if(any(abs(nnC - nnI) > eps)) stop("Algorithms for nndist() do not agree") nn3C <- nndist(X, k=3, method="C") nn3I <- nndist(X, k=3, method="interpreted") if(any(abs(nn3C - nn3I) > eps)) stop("Algorithms for nndist(k=3) do not agree") # nnwhich.ppp nwC <- nnwhich(X, method="C") nwI <- nnwhich(X, method="interpreted") if(any(nwC != nwI)) stop("Algorithms for nnwhich() do not agree") nw3C <- nnwhich(X, k=3, method="C") nw3I <- nnwhich(X, k=3, method="interpreted") if(any(nw3C != nw3I)) stop("Algorithms for nnwhich(k=3) do not agree") }) # ppmBadData.R # $Revision: 1.4 $ $Date: 2011/12/05 07:29:16 $ # Testing robustness of ppm and support functions # when data are rubbish require(spatstat) local({ # --------------------------------------------------- # from Rolf: very large proportion of data is NA SEED <- 42 K <- 101 A <- 500 X <- seq(0, A, length=K) G <- expand.grid(x=X, y=X) FOO <- function(x,y) { sin(x)^2 + cos(y)^2 } M1 <- im(matrix(FOO(G$x, G$y), K, K), xcol=X, yrow=X) M <- im(matrix(FOO(G$x, G$y), K, K)) BAR <- function(x) { exp(-6.618913 + 5.855337 * x - 8.432483 * x^2) } V <- im(BAR(M$v), xcol=X, yrow=X) # V <- eval.im(exp(-6.618913 + 5.855337 * M - 8.432483 * M^2)) set.seed(SEED) Y <- rpoispp(V) fY <- ppm(Y, ~cv + I(cv^2), covariates=list(cv=M), correction="translate") diagnose.ppm(fY) lurking(fY, covariate=as.im(function(x,y){x}, square(A)), type="raw") }) # -------------------------------------------------------- # from Andrew Bevan: numerical overflow, ill-conditioned Fisher information local({ SEED <- 42 nongranite<- owin(poly = list(x = c(0, 8500, 7000, 6400, 6400, 6700, 7000, 7200, 7300, 8000, 8100, 8800, 9500, 10000, 10000, 0), y = c(0, 0, 2000, 3800, 4000, 5000, 6500, 7400, 7500, 8000, 8100, 9000, 9500, 9600, 10000, 10000))) #Trend on raster grid rain <- as.im(X=function(x,y) { x^2 + y^2 }, W=nongranite, dimyx=100) #Generate a point pattern via a Lennard-Jones process set.seed(SEED) mod4<- rmhmodel(cif="lennard", par=list(beta=1, sigma=250, epsilon=2.2), trend=rain, w=nongranite) ljtr<- rmh(mod4, start=list(n.start=80), control=list(p=1, nrep=1e5)) #Fit a point process model to the pattern with rain as a covariate # NOTE INCORRECT TREND FORMULA ljtrmod <- ppm(ljtr, trend= ~ Z, interaction=NULL, covariates=list(Z=rain)) ss <- summary(ljtrmod) }) # # tests/ppmgam.R # # Test ppm with use.gam=TRUE # # $Revision: 1.2 $ $Date: 2013/03/01 08:46:34 $ # require(spatstat) local({ fit <- ppm(nztrees, ~s(x,y), use.gam=TRUE) mm <- model.matrix(fit) mf <- model.frame(fit) v <- vcov(fit) }) # ppmmarkorder.R # $Revision: 1.2 $ $Date: 2011/12/05 07:29:16 $ # Test that predict.ppm, plot.ppm and plot.fitin # tolerate marks with levels that are not in alpha order # require(spatstat) local({ X <- amacrine levels(marks(X)) <- c("ZZZ", "AAA") fit <- ppm(X, ~marks, MultiStrauss(c("ZZZ","AAA"), matrix(0.06, 2, 2))) aa <- predict(fit, type="trend") bb <- predict(fit, type="cif") plot(fit) plot(fitin(fit)) }) # # tests/ppmscope.R # # Test things that might corrupt the internal format of ppm objects # # $Revision: 1.3 $ $Date: 2011/12/05 07:29:16 $ # # (1) Scoping problem that can arise when ppm splits the data require(spatstat) local({ fit <- ppm(bei, ~elev, covariates=bei.extra) mm <- model.matrix(fit) # (2) Fast update mechanism fit1 <- ppm(cells, ~x+y, Strauss(0.07)) fit2 <- update(fit1, ~y) fit3 <- update(fit2, ~x) }) # # tests/rmhAux.R # # $Revision: 1.1 $ $Date: 2013/02/18 10:41:27 $ # # For interactions which maintain 'auxiliary data', # verify that the auxiliary data are correctly updated. # # To do this we run rmh with nsave=1 so that the point pattern state # is saved after every iteration, then the algorithm is restarted, # and the auxiliary data are re-initialised. The final state must agree with # the result of simulation without saving. # ---------------------------------------------------- require(spatstat) local({ # Geyer: mod <- list(cif="geyer", par=list(beta=1.25,gamma=1.6,r=0.2,sat=4.5), w=square(10)) set.seed(42) X.nosave <- rmh(model=mod, start=list(n.start=50), control=list(nrep=1e3, periodic=FALSE, expand=1)) set.seed(42) X.save <- rmh(model=mod, start=list(n.start=50), control=list(nrep=1e3, periodic=FALSE, expand=1, nburn=0, nsave=1)) stopifnot(npoints(X.save) == npoints(X.nosave)) stopifnot(max(nncross(X.save, X.nosave)$dist) == 0) stopifnot(max(nncross(X.nosave, X.save)$dist) == 0) }) # Test examples for rmh.default # run to reasonable length # and with tests for validity added # ---------------------------------------------------- require(spatstat) local({ if(!exists("nr")) nr <- 5e3 spatstat.options(expand=1.1) # Strauss process. mod01 <- list(cif="strauss",par=list(beta=2,gamma=0.2,r=0.7), w=c(0,10,0,10)) X1.strauss <- rmh(model=mod01,start=list(n.start=80), control=list(nrep=nr)) # Strauss process, conditioning on n = 80: X2.strauss <- rmh(model=mod01,start=list(n.start=80), control=list(p=1,nrep=nr)) stopifnot(X2.strauss$n == 80) # test tracking mechanism X1.strauss <- rmh(model=mod01,start=list(n.start=80), control=list(nrep=nr), track=TRUE) X2.strauss <- rmh(model=mod01,start=list(n.start=80), control=list(p=1,nrep=nr), track=TRUE) # Hard core process: mod02 <- list(cif="hardcore",par=list(beta=2,hc=0.7),w=c(0,10,0,10)) X3.hardcore <- rmh(model=mod02,start=list(n.start=60), control=list(nrep=nr)) # Strauss process equal to pure hardcore: mod02 <- list(cif="strauss",par=list(beta=2,gamma=0,r=0.7),w=c(0,10,0,10)) X3.strauss <- rmh(model=mod02,start=list(n.start=60), control=list(nrep=nr)) # Strauss process in a polygonal window. x <- c(0.55,0.68,0.75,0.58,0.39,0.37,0.19,0.26,0.42) y <- c(0.20,0.27,0.68,0.99,0.80,0.61,0.45,0.28,0.33) mod03 <- list(cif="strauss",par=list(beta=2000,gamma=0.6,r=0.07), w=owin(poly=list(x=x,y=y))) X4.strauss <- rmh(model=mod03,start=list(n.start=90), control=list(nrep=nr)) # Strauss process in a polygonal window, conditioning on n = 42. X5.strauss <- rmh(model=mod03,start=list(n.start=42), control=list(p=1,nrep=nr)) stopifnot(X5.strauss$n == 42) # Strauss process, starting off from X4.strauss, but with the # polygonal window replace by a rectangular one. At the end, # the generated pattern is clipped to the original polygonal window. xxx <- X4.strauss xxx$window <- as.owin(c(0,1,0,1)) X6.strauss <- rmh(model=mod03,start=list(x.start=xxx), control=list(nrep=nr)) # Strauss with hardcore: mod04 <- list(cif="straush",par=list(beta=2,gamma=0.2,r=0.7,hc=0.3), w=c(0,10,0,10)) X1.straush <- rmh(model=mod04,start=list(n.start=70), control=list(nrep=nr)) # Another Strauss with hardcore (with a perhaps surprising result): mod05 <- list(cif="straush",par=list(beta=80,gamma=0.36,r=45,hc=2.5), w=c(0,250,0,250)) X2.straush <- rmh(model=mod05,start=list(n.start=250), control=list(nrep=nr)) # Pure hardcore (identical to X3.strauss). mod06 <- list(cif="straush",par=list(beta=2,gamma=1,r=1,hc=0.7), w=c(0,10,0,10)) X3.straush <- rmh(model=mod06,start=list(n.start=60), control=list(nrep=nr)) # Area-interaction, inhibitory mod.area <- list(cif="areaint",par=list(beta=2,eta=0.5,r=0.5), w=square(10)) X.area <- rmh(model=mod.area,start=list(n.start=60), control=list(nrep=nr)) # Area-interaction, clustered mod.area2 <- list(cif="areaint",par=list(beta=2,eta=1.5,r=0.5), w=square(10)) X.area2 <- rmh(model=mod.area2,start=list(n.start=60), control=list(nrep=nr)) # Area-interaction close to hard core set.seed(42) mod.area0 <- list(cif="areaint",par=list(beta=2,eta=1e-300,r=0.35), w=square(10)) X.area0 <- rmh(model=mod.area0,start=list(x.start=X3.hardcore), control=list(nrep=nr)) stopifnot(nndist(X.area0) > 0.6) # Soft core: w <- c(0,10,0,10) mod07 <- list(cif="sftcr",par=list(beta=0.8,sigma=0.1,kappa=0.5), w=c(0,10,0,10)) X.sftcr <- rmh(model=mod07,start=list(n.start=70), control=list(nrep=nr)) # Diggle, Gates, and Stibbard: mod12 <- list(cif="dgs",par=list(beta=3600,rho=0.08),w=c(0,1,0,1)) X.dgs <- rmh(model=mod12,start=list(n.start=300), control=list(nrep=nr)) # Diggle-Gratton: mod13 <- list(cif="diggra", par=list(beta=1800,kappa=3,delta=0.02,rho=0.04), w=square(1)) X.diggra <- rmh(model=mod13,start=list(n.start=300), control=list(nrep=nr)) # Geyer: mod14 <- list(cif="geyer",par=list(beta=1.25,gamma=1.6,r=0.2,sat=4.5), w=c(0,10,0,10)) X1.geyer <- rmh(model=mod14,start=list(n.start=200), control=list(nrep=nr)) # Geyer; same as a Strauss process with parameters # (beta=2.25,gamma=0.16,r=0.7): mod15 <- list(cif="geyer",par=list(beta=2.25,gamma=0.4,r=0.7,sat=10000), w=c(0,10,0,10)) X2.geyer <- rmh(model=mod15,start=list(n.start=200), control=list(nrep=nr)) mod16 <- list(cif="geyer",par=list(beta=8.1,gamma=2.2,r=0.08,sat=3)) data(redwood) X3.geyer <- rmh(model=mod16,start=list(x.start=redwood), control=list(periodic=TRUE,nrep=nr)) # Geyer, starting from the redwood data set, simulating # on a torus, and conditioning on n: X4.geyer <- rmh(model=mod16,start=list(x.start=redwood), control=list(p=1,periodic=TRUE,nrep=nr)) # Lookup (interaction function h_2 from page 76, Diggle (2003)): r <- seq(from=0,to=0.2,length=101)[-1] # Drop 0. h <- 20*(r-0.05) h[r<0.05] <- 0 h[r>0.10] <- 1 mod17 <- list(cif="lookup",par=list(beta=4000,h=h,r=r),w=c(0,1,0,1)) X.lookup <- rmh(model=mod17,start=list(n.start=100), control=list(nrep=nr)) # Strauss with trend tr <- function(x,y){x <- x/250; y <- y/250; exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6) } beta <- 0.3 gmma <- 0.5 r <- 45 tr3 <- function(x,y){x <- x/250; y <- y/250; exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6) } # log quadratic trend mod17 <- list(cif="strauss",par=list(beta=beta,gamma=gmma,r=r),w=c(0,250,0,250), trend=tr3) X1.strauss.trend <- rmh(model=mod17,start=list(n.start=90), control=list(nrep=nr)) }) # Things which should cause an error require(spatstat) local({ if(!exists("nv")) nv <- 0 if(!exists("nr")) nr <- 1e3 # Strauss with zero intensity and p = 1 mod0S <- list(cif="strauss",par=list(beta=0,gamma=0.6,r=0.7), w = square(3)) out <- try(X0S <- rmh(model=mod0S,start=list(n.start=80), control=list(p=1,nrep=nr,nverb=nv),verbose=FALSE)) if(!inherits(out, "try-error")) stop("Error not trapped (Strauss with zero intensity and p = 1) in tests/rmhErrors.R") }) # # tests/rmhExpand.R # # test decisions about expansion of simulation window # # $Revision: 1.2 $ $Date: 2011/12/05 07:29:16 $ # require(spatstat) local({ fit <- ppm(cells, ~x) # check rmhmodel.ppm mod <- rmhmodel(fit) wsim <- as.rectangle(mod$trend) if(!identical(wsim, as.owin(cells))) stop("Expansion occurred improperly in rmhmodel.ppm") }) # # tests of rmh, running multitype point processes # require(spatstat) local({ if(!exists("nr")) nr <- 5e3 if(!exists("nv")) nv <- 0 spatstat.options(expand=1.1) # Multitype Poisson modp2 <- list(cif="poisson", par=list(beta=2), types=letters[1:3], w = square(10)) Xp2 <- rmh(modp2, start=list(n.start=0), control=list(p=1)) # Multitype Strauss: beta <- c(0.027,0.008) gmma <- matrix(c(0.43,0.98,0.98,0.36),2,2) r <- matrix(c(45,45,45,45),2,2) mod08 <- list(cif="straussm",par=list(beta=beta,gamma=gmma,radii=r), w=c(0,250,0,250)) X1.straussm <- rmh(model=mod08,start=list(n.start=80), control=list(ptypes=c(0.75,0.25),nrep=nr,nverb=nv)) # Multitype Strauss conditioning upon the total number # of points being 80: X2.straussm <- rmh(model=mod08,start=list(n.start=80), control=list(p=1,ptypes=c(0.75,0.25),nrep=nr, nverb=nv)) stopifnot(X2.straussm$n == 80) # Conditioning upon the number of points of type 1 being 60 # and the number of points of type 2 being 20: X3.straussm <- rmh(model=mod08,start=list(n.start=c(60,20)), control=list(fixall=TRUE,p=1,ptypes=c(0.75,0.25), nrep=nr,nverb=nv)) stopifnot(all(table(X3.straussm$marks) == c(60,20))) # Multitype Strauss hardcore: rhc <- matrix(c(9.1,5.0,5.0,2.5),2,2) mod09 <- list(cif="straushm",par=list(beta=beta,gamma=gmma, iradii=r,hradii=rhc),w=c(0,250,0,250)) X.straushm <- rmh(model=mod09,start=list(n.start=80), control=list(ptypes=c(0.75,0.25),nrep=nr,nverb=nv)) # Multitype Strauss hardcore with trends for each type: beta <- c(0.27,0.08) tr3 <- function(x,y){x <- x/250; y <- y/250; exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6) } # log quadratic trend tr4 <- function(x,y){x <- x/250; y <- y/250; exp(-0.6*x+0.5*y)} # log linear trend mod10 <- list(cif="straushm",par=list(beta=beta,gamma=gmma, iradii=r,hradii=rhc),w=c(0,250,0,250), trend=list(tr3,tr4)) X1.straushm.trend <- rmh(model=mod10,start=list(n.start=350), control=list(ptypes=c(0.75,0.25), nrep=nr,nverb=nv)) # Multitype Strauss hardcore with trends for each type, given as images: bigwin <- square(250) i1 <- as.im(tr3, bigwin) i2 <- as.im(tr4, bigwin) mod11 <- list(cif="straushm",par=list(beta=beta,gamma=gmma, iradii=r,hradii=rhc),w=bigwin, trend=list(i1,i2)) X2.straushm.trend <- rmh(model=mod11,start=list(n.start=350), control=list(ptypes=c(0.75,0.25),expand=1, nrep=nr,nverb=nv)) ####################################################################### ############ checks on distribution of output ####################### ####################################################################### checkp <- function(p, context, testname, failmessage, pcrit=0.01) { if(missing(failmessage)) failmessage <- paste("output failed", testname) if(p < pcrit) warning(paste(context, ",", failmessage), call.=FALSE) cat(paste("\n", context, ",", testname, "has p-value", signif(p,4), "\n")) } # Multitype Strauss code; output is multitype Poisson beta <- 100 * c(1,1) ri <- matrix(0.07, 2, 2) gmma <- matrix(1, 2, 2) # no interaction tr1 <- function(x,y){ rep(1, length(x)) } tr2 <- function(x,y){ rep(2, length(x)) } mod <- rmhmodel(cif="straussm", par=list(beta=beta,gamma=gmma,radii=ri), w=owin(), trend=list(tr1,tr2)) X <- rmh(mod, start=list(n.start=0), control=list(nrep=1e6)) # The model is Poisson with intensity 100 for type 1 and 200 for type 2. # Total number of points is Poisson (300) # Marks are i.i.d. with P(type 1) = 1/3, P(type 2) = 2/3. # Test whether the total intensity looks right # p <- ppois(X$n, 300) p.val <- 2 * min(p, 1-p) checkp(p.val, "In multitype Poisson simulation", "test whether total number of points has required mean value") # Test whether the mark distribution looks right ta <- table(X$marks) cat("Frequencies of marks:") print(ta) checkp(chisq.test(ta, p = c(1,2)/3)$p.value, "In multitype Poisson simulation", "chi-squared goodness-of-fit test for mark distribution (1/3, 2/3)") ##### #### multitype Strauss code; fixall=TRUE; #### output is multinomial process with nonuniform locations #### the.context <- "In nonuniform multinomial simulation" beta <- 100 * c(1,1) ri <- matrix(0.07, 2, 2) gmma <- matrix(1, 2, 2) # no interaction tr1 <- function(x,y){ ifelse(x < 0.5, 0, 2) } tr2 <- function(x,y){ ifelse(y < 0.5, 1, 3) } # cdf of these distributions Fx1 <- function(x) { ifelse(x < 0.5, 0, ifelse(x < 1, 2 * x - 1, 1)) } Fy2 <- function(y) { ifelse(y < 0, 0, ifelse(y < 0.5, y/2, ifelse(y < 1, (1/2 + 3 * (y-1/2))/2, 1))) } mod <- rmhmodel(cif="straussm", par=list(beta=beta,gamma=gmma,radii=ri), w=owin(), trend=list(tr1,tr2)) X <- rmh(mod, start=list(n.start=c(50,50)), control=list(nrep=1e6, expand=1, p=1, fixall=TRUE)) # The model is Poisson # Mean number of type 1 points = 100 # Mean number of type 2 points = 200 # Total intensity = 300 # Marks are i.i.d. with P(type 1) = 1/3, P(type 2) = 2/3 # Test whether the coordinates look OK Y <- split(X) X1 <- Y[[names(Y)[1]]] X2 <- Y[[names(Y)[2]]] checkp(ks.test(X1$y, "punif")$p.value, the.context, "Kolmogorov-Smirnov test of uniformity of y coordinates of type 1 points") if(any(X1$x < 0.5)) { warning(paste(the.context, ",", "x-coordinates of type 1 points are IMPOSSIBLE"), call.=FALSE) } else { checkp(ks.test(Fx1(X1$x), "punif")$p.value, the.context, "Kolmogorov-Smirnov test of uniformity of transformed x coordinates of type 1 points") } checkp(ks.test(X2$x, "punif")$p.value, the.context, "Kolmogorov-Smirnov test of uniformity of x coordinates of type 2 points") checkp(ks.test(Fy2(X2$y), "punif")$p.value, the.context, "Kolmogorov-Smirnov test of uniformity of transformed y coordinates of type 2 points") }) # # tests/rmhTrend.R # # Problems with trend images (rmhmodel.ppm or rmhEngine) # require(spatstat) local({ set.seed(42) # Bug folder 37 of 8 feb 2011 # rmhmodel.ppm -> predict.ppm # + rmhResolveTypes -> is.subset.owin data(demopat) Z <- rescale(demopat, 7000) X <- unmark(Z) X1 <- split(Z)[[1]] Int <- density(X,dimyx=200) Lint <- eval.im(log(npoints(X1)*Int/npoints(X))) M <- as.owin(Int) MR <- intersect.owin(M,scalardilate(M,0.5,origin="midpoint")) X1 <- X1[MR] Fut <- ppm(X1,~offset(Lint),covariates=list(Lint=Lint), inter=BadGey(r=c(0.03,0.05),sat=3)) Y <- rmh(Fut,control=list(expand=M,nrep=1e3), verbose=FALSE) }) # strange boundary cases require(spatstat) local({ if(!exists("nv")) nv <- 0 if(!exists("nr")) nr <- 5e3 # Poisson process cat("Poisson\n") modP <- list(cif="poisson",par=list(beta=10), w = square(3)) XP <- rmh(model = modP, start = list(n.start=25), control=list(nrep=nr,nverb=nv)) # Poisson process case of Strauss cat("\nPoisson case of Strauss\n") modPS <- list(cif="strauss",par=list(beta=10,gamma=1,r=0.7), w = square(3)) XPS <- rmh(model=modPS, start=list(n.start=25), control=list(nrep=nr,nverb=nv)) # Strauss with zero intensity cat("\nStrauss with zero intensity\n") mod0S <- list(cif="strauss",par=list(beta=0,gamma=0.6,r=0.7), w = square(3)) X0S <- rmh(model=mod0S,start=list(n.start=80), control=list(nrep=nr,nverb=nv)) stopifnot(X0S$n == 0) # Poisson with zero intensity cat("\nPoisson with zero intensity\n") mod0P <- list(cif="poisson",par=list(beta=0), w = square(3)) X0P <- rmh(model = mod0P, start = list(n.start=25), control=list(nrep=nr,nverb=nv)) # Poisson conditioned on zero points cat("\nPoisson conditioned on zero points\n") modp <- list(cif="poisson", par=list(beta=2), w = square(10)) Xp <- rmh(modp, start=list(n.start=0), control=list(p=1, nrep=nr)) stopifnot(Xp$n == 0) # Multitype Poisson conditioned on zero points cat("\nMultitype Poisson conditioned on zero points\n") modp2 <- list(cif="poisson", par=list(beta=2), types=letters[1:3], w = square(10)) Xp2 <- rmh(modp2, start=list(n.start=0), control=list(p=1, nrep=nr)) stopifnot(is.marked(Xp2)) stopifnot(Xp2$n == 0) # Multitype Poisson conditioned on zero points of each type cat("\nMultitype Poisson conditioned on zero points of each type\n") Xp2fix <- rmh(modp2, start=list(n.start=c(0,0,0)), control=list(p=1, fixall=TRUE, nrep=nr)) stopifnot(is.marked(Xp2fix)) stopifnot(Xp2fix$n == 0) }) # # tests of rmhmodel.ppm # require(spatstat) local({ f <- ppm(cells) m <- rmhmodel(f) f <- ppm(cells, ~x) m <- rmhmodel(f) f <- ppm(cells, ~1, Strauss(0.1)) m <- rmhmodel(f) f <- ppm(cells, ~1, StraussHard(r=0.1,hc=0.05)) m <- rmhmodel(f) f <- ppm(cells, ~1, Hardcore(0.07)) m <- rmhmodel(f) f <- ppm(cells, ~1, DiggleGratton(0.05,0.1)) m <- rmhmodel(f) f <- ppm(cells, ~1, Softcore(0.5), correction="isotropic") m <- rmhmodel(f) f <- ppm(cells, ~1, Geyer(0.07,2)) m <- rmhmodel(f) f <- ppm(cells, ~1, BadGey(c(0.07,0.1,0.13),2)) m <- rmhmodel(f) f <- ppm(cells, ~1, PairPiece(r = c(0.05, 0.1, 0.2))) m <- rmhmodel(f) f <- ppm(cells, ~1, AreaInter(r=0.06)) m <- rmhmodel(f) # multitype r <- matrix(0.07, 2, 2) f <- ppm(amacrine, ~1, MultiStrauss(c("off","on"),r)) m <- rmhmodel(f) h <- matrix(min(nndist(amacrine))/2, 2, 2) f <- ppm(amacrine, ~1, MultiStraussHard(c("off","on"),r, h)) m <- rmhmodel(f) diag(r) <- NA diag(h) <- NA f <- ppm(amacrine, ~1, MultiStrauss(c("off","on"),r)) m <- rmhmodel(f) f <- ppm(amacrine, ~1, MultiStraussHard(c("off","on"),r, h)) m <- rmhmodel(f) # multitype data, interaction not dependent on type f <- ppm(amacrine, ~marks, Strauss(0.05)) m <- rmhmodel(f) # trends f <- ppm(cells, ~x, Strauss(0.1)) m <- rmhmodel(f) f <- ppm(cells, ~y, StraussHard(r=0.1,hc=0.05)) m <- rmhmodel(f) f <- ppm(cells, ~x+y, Hardcore(0.07)) m <- rmhmodel(f) f <- ppm(cells, ~polynom(x,y,2), Softcore(0.5), correction="isotropic") m <- rmhmodel(f) # covariates Z <- as.im(function(x,y){ x^2+y^2 }, as.owin(cells)) f <- ppm(cells, ~z, covariates=list(z=Z)) m <- rmhmodel(f) m <- rmhmodel(f, control=list(p=1)) Zim <- as.im(Z, as.owin(cells)) f <- ppm(amacrine, ~z, covariates=list(z=Zim)) m <- rmhmodel(f) Z <- as.im(function(x,y){ x^2+y }, as.owin(amacrine)) f <- ppm(amacrine, ~z + marks, covariates=list(z=Z)) m <- rmhmodel(f) m <- rmhmodel(f, control=list(p=1)) m <- rmhmodel(f, control=list(p=1,fixall=TRUE)) Zim <- as.im(Z, as.owin(amacrine)) f <- ppm(amacrine, ~z + marks, covariates=list(z=Zim)) m <- rmhmodel(f) }) # # tests/rmhmodelHybrids.R # # Test that rmhmodel.ppm and rmhmodel.default # work on Hybrid interaction models # # $Revision: 1.3 $ $Date: 2013/01/29 02:12:13 $ # require(spatstat) local({ # ......... rmhmodel.ppm ....................... fit1 <- ppm(redwood, ~1, Hybrid(A=Strauss(0.02), B=Geyer(0.1, 2), C=Geyer(0.15, 1))) m1 <- rmhmodel(fit1) m1 reach(m1) # Test of handling 'IsOffset' data(cells) fit2 <- ppm(cells, ~1, Hybrid(H=Hardcore(0.05), G=Geyer(0.15, 2))) rmhmodel(fit2) # Test of handling Poisson components fit3 <- ppm(cells, ~1, Hybrid(P=Poisson(), S=Strauss(0.05))) X3 <- rmh(fit3, control=list(nrep=1e3,expand=1), verbose=FALSE) # ............ rmhmodel.default ............................ modH <- list(cif=c("strauss","geyer"), par=list(list(beta=50,gamma=0.5, r=0.1), list(beta=1, gamma=0.7, r=0.2, sat=2)), w = square(1)) rmodH <- rmhmodel(modH) rmodH reach(rmodH) # test handling of Poisson components modHP <- list(cif=c("poisson","strauss"), par=list(list(beta=5), list(beta=10,gamma=0.5, r=0.1)), w = square(1)) rmodHP <- rmhmodel(modHP) rmodHP reach(rmodHP) modPP <- list(cif=c("poisson","poisson"), par=list(list(beta=5), list(beta=10)), w = square(1)) rmodPP <- rmhmodel(modPP) rmodPP reach(rmodPP) }) # # tests/rmh.ppm.R # # $Revision: 1.1 $ $Date: 2012/10/14 07:24:21 $ # # Examples removed from rmh.ppm.Rd # stripped down to minimal tests of validity # require(spatstat) local({ op <- spatstat.options() spatstat.options(rmh.nrep=10, npixel=10, ndummy.min=10) spatstat.options(project.fast=TRUE) Nrep <- 10 X <- swedishpines # Poisson process fit <- ppm(X, ~1, Poisson()) Xsim <- rmh(fit) # Strauss process fit <- ppm(X, ~1, Strauss(r=7)) Xsim <- rmh(fit) # Strauss process simulated on a larger window # then clipped to original window Xsim <- rmh(fit, control=list(nrep=Nrep, expand=1.1, periodic=TRUE)) # Strauss - hard core process # fit <- ppm(X, ~1, StraussHard(r=7,hc=2)) # Xsim <- rmh(fit, start=list(n.start=X$n)) # Geyer saturation process # fit <- ppm(X, ~1, Geyer(r=7,sat=2)) # Xsim <- rmh(fit, start=list(n.start=X$n)) # Area-interaction process fit <- ppm(X, ~1, AreaInter(r=7)) Xsim <- rmh(fit, start=list(n.start=X$n)) # soft core interaction process # X <- quadscheme(X, nd=50) # fit <- ppm(X, ~1, Softcore(kappa=0.1), correction="isotropic") # Xsim <- rmh(fit, start=list(n.start=X$n)) # Diggle-Gratton pairwise interaction model # fit <- ppm(cells, ~1, DiggleGratton(0.05, 0.1)) # Xsim <- rmh(fit, start=list(n.start=cells$n)) # plot(Xsim, main="simulation from fitted Diggle-Gratton model") X <- rSSI(0.05, 100) # piecewise-constant pairwise interaction function fit <- ppm(X, ~1, PairPiece(seq(0.02, 0.1, by=0.01))) Xsim <- rmh(fit) # marked point pattern Y <- amacrine # marked Poisson models fit <- ppm(Y) Ysim <- rmh(fit) fit <- ppm(Y,~marks) Ysim <- rmh(fit) fit <- ppm(Y,~x) Ysim <- rmh(fit) # fit <- ppm(Y,~polynom(x,2)) # Ysim <- rmh(fit) fit <- ppm(Y,~marks+x) Ysim <- rmh(fit) # fit <- ppm(Y,~marks+polynom(x,2)) # Ysim <- rmh(fit) # multitype Strauss models MS <- MultiStrauss(types = levels(Y$marks), radii=matrix(0.07, ncol=2, nrow=2)) # fit <- ppm(Y,~marks*polynom(x,2), MS) fit <- ppm(Y,~marks*x, MS) Ysim <- rmh(fit) spatstat.options(op) }) # fvproblems.R require(spatstat) # This appears in the workshop notes # Problem detected by Martin Bratschi local({ Jdif <- function(X, ..., i) { Jidot <- Jdot(X, ..., i=i) J <- Jest(X, ...) dif <- eval.fv(Jidot - J) return(dif) } Z <- Jdif(amacrine, i="on") }) # # Test mathlegend code # local({ K <- Kest(cells) plot(K) plot(K, . ~ r) plot(K, . - theo ~ r) plot(K, sqrt(./pi) ~ r) plot(K, cbind(iso, theo) ~ r) plot(K, cbind(iso, theo) - theo ~ r) plot(K, sqrt(cbind(iso, theo)/pi) ~ r) plot(K, cbind(iso/2, -theo) ~ r) plot(K, cbind(iso/2, trans/2) - theo ~ r) # test expansion of .x and .y plot(K, . ~ .x) plot(K, . - theo ~ .x) plot(K, .y - theo ~ .x) plot(K, sqrt(.y) - sqrt(theo) ~ .x) # problems with parsing weird strings in levels(marks(X)) # noted by Ulf Mehlig levels(marks(amacrine)) <- c("Nastricreechia krorluppia", "Homo habilis") plot(Kcross(amacrine)) plot(alltypes(amacrine, "K")) plot(alltypes(amacrine, "J")) plot(alltypes(amacrine, pcfcross)) }) # test of case where mark levels contain illegal characters require(spatstat) local({ hyphenated <- c("a", "not-a") spaced <- c("U", "non U") suffixed <- c("a+", "a*") charred <- c("+", "*") irad <- matrix(0.1, 2,2) hrad <- matrix(0.005, 2, 2) tryit <- function(types, X, irad, hrad) { levels(marks(X)) <- types fit <- ppm(X, ~marks + polynom(x,y,2), MultiStraussHard(types=types,iradii=irad,hradii=hrad)) print(fit) print(coef(fit)) val <- fitted(fit) pred <- predict(fit) return(invisible(NULL)) } tryit(hyphenated, amacrine, irad, hrad) tryit(spaced, amacrine, irad, hrad) tryit(suffixed, amacrine, irad, hrad) tryit(charred, amacrine, irad, hrad) }) # # test cases where there are no (rows or columns of) marks # require(spatstat) local({ n <- npoints(cells) df <- data.frame(x=1:n, y=factor(sample(letters, n, replace=TRUE))) nocolumns <- c(FALSE, FALSE) norows <- rep(FALSE, n) X <- cells marks(X) <- df marks(X) <- df[,1] marks(X) <- df[,nocolumns] Z <- Y <- X[integer(0)] marks(Y) <- df[norows,] stopifnot(is.marked(Y)) marks(Z) <- df[norows,nocolumns] stopifnot(!is.marked(Z)) }) # checks validity of fast C implementation of Geyer interaction require(spatstat) local({ X <- redwood Q <- quadscheme(X) U <- union.quad(Q) EP <- equalpairs.quad(Q) G <- Geyer(0.11, 2) # The value r=0.11 is chosen to avoid hardware numerical effects (gcc bug 323). # It avoids being close any value of pairdist(redwood). # The nearest such values are 0.1077.. and 0.1131.. # By contrast if r = 0.1 there are values differing from 0.1 by 3e-17 a <- pairsat.family$eval(X,U,EP,G$pot,G$par,"border") b <- G$fasteval(X,U,EP,G$pot,G$par,"border") if(!all(a==b)) stop("Results of Geyer()$fasteval and pairsat.family$eval do not match") # ... # and again for a non-integer value of 'sat' # (spotted by Thordis Linda Thorarinsdottir) G <- Geyer(0.11, 2.5) a <- pairsat.family$eval(X,U,EP,G$pot,G$par,"border") b <- G$fasteval(X,U,EP,G$pot,G$par,"border") if(!all(a==b)) stop("Results of Geyer()$fasteval and pairsat.family$eval do not match when sat is not an integer") }) require(spatstat) local({ co <- as.ppp(corners(letterR), letterR, check=FALSE) co[letterR] }) # tests/segments.R # $Revision: 1.7 $ $Date: 2011/12/05 07:29:16 $ require(spatstat) local({ # pointed out by Jeff Laake W <- owin() X <- psp(x0=.25,x1=.25,y0=0,y1=1,window=W) X[W] # migrated from 'lpp' X <- psp(runif(10),runif(10),runif(10),runif(10), window=owin()) Z <- as.mask.psp(X) Z <- pixellate(X) # test of distppll pointed out by Ang Qi Wei p <- matrix(c(1.5, 0), 1, 2) l <- matrix(c(0,0,1,0,1,0,2,0), 2, 4, byrow=T) a <- distppll(p, l, mintype=2, method="interpreted") b <- distppll(p, l, mintype=2, method="Fortran") d <- distppll(p, l, mintype=2, method="C") if(a$min.which != b$min.which) stop("conflict between Fortran and interpreted code in distppll") if(a$min.which != d$min.which) stop("conflict between C and interpreted code in distppll") # tests of pixellate.psp -> seg2pixL ns <- 50 out <- numeric(ns) for(i in 1:ns) { X <- psp(runif(1), runif(1), runif(1), runif(1), window=owin()) len <- lengths.psp(X) dlen <- sum(pixellate(X)$v) out[i] <- if(len > 1e-7) dlen/len else 1 } if(diff(range(out)) > 0.01) stop(paste( "pixellate.psp test 1: relative error [", paste(diff(range(out)), collapse=", "), "]")) # Michael Sumner's test examples set.seed(33) n <- 2001 co <- cbind(runif(n), runif(n)) ow <- owin() X <- psp(co[-n,1], co[-n,2], co[-1,1], co[-1,2], window=ow) s1 <- sum(pixellate(X)) s2 <- sum(lengths.psp(X)) if(abs(s1 - s2)/s2 > 0.01) { stop(paste("pixellate.psp test 2:", "sum(pixellate(X)) = ", s1, "!=", s2, "= sum(lengths.psp(X))")) } wts <- 1/(lengths.psp(X) * X$n) s1 <- sum(pixellate(X, weights=wts)) if(abs(s1-1) > 0.01) { stop(paste("pixellate.psp test 3:", "sum(pixellate(X, weights))=", s1, " (should be 1)")) } X <- psp(0, 0, 0.01, 0.001, window=owin()) s1 <- sum(pixellate(X)) s2 <- sum(lengths.psp(X)) if(abs(s1 - s2)/s2 > 0.01) { stop(paste("pixellate.psp test 4:", "sum(pixellate(X)) = ", s1, "!=", s2, "= sum(lengths.psp(X))")) } X <- psp(0, 0, 0.001, 0.001, window=owin()) s1 <- sum(pixellate(X)) s2 <- sum(lengths.psp(X)) if(abs(s1 - s2)/s2 > 0.01) { stop(paste("pixellate.psp test 5:", "sum(pixellate(X)) = ", s1, "!=", s2, "= sum(lengths.psp(X))")) } }) # test for step() operation # require(spatstat) local({ Z <- as.im(function(x,y){ x^3 - y^2 }, nztrees$window) fitP <- ppm(nztrees, ~x+y+Z, covariates=list(Z=Z)) step(fitP) fitS <- update(fitP, Strauss(7)) step(fitS) fitM <- ppm(amacrine, ~ marks*(x+y), MultiStrauss(types=levels(marks(amacrine)), radii=matrix(0.04, 2, 2))) step(fitM) }) require(spatstat) source("badwindow.R") owinpolycheck(W,verbose=FALSE) # # tests/hobjects.R # # Validity of methods for ppm(... method="ho") # require(spatstat) local({ fit <- ppm(cells, ~1, Strauss(0.1), method="ho", nsim=10) fitx <- ppm(cells, ~offset(x), Strauss(0.1), method="ho", nsim=10) a <- AIC(fit) ax <- AIC(fitx) f <- fitted(fit) fx <- fitted(fitx) p <- predict(fit) px <- predict(fitx) }) # check kstest with strange data require(spatstat) local({ # Marked point patterns with some marks not represented AC <- split(ants, un=FALSE)$Cataglyphis AM <- split(ants, un=FALSE)$Messor DM <- distmap(AM) # should produce a warning, rather than a crash: kstest(AC, DM) # should be OK: kstest(unmark(AC), DM) # linear networks X <- runiflpp(20, simplenet) fit <- lppm(X, ~1) kstest(fit, "y") }) # # tests/kppm.R # # $Revision: 1.6 $ $Date: 2012/04/08 03:22:20 $ # # Test functionality of kppm that depends on RandomFields # require(spatstat) local({ if(require(RandomFields) && RandomFieldsSafe()) { fit0 <- kppm(redwood, ~1, "LGCP") simulate(fit0) fit <- kppm(redwood, ~x, "LGCP", covmodel=list(model="matern", nu=0.3), control=list(maxit=5)) simulate(fit) # ... and Abdollah's code fit <- kppm(redwood, ~x, cluster="Cauchy", statistic="K") simulate(fit) } }) # temporary test file for localpcfmatrix require(spatstat) local({ a <- localpcfmatrix(redwood) a plot(a) a[, 3:5] }) # check for various bugs related to factor conversions require(spatstat) local({ # make a factor image m <- factor(rep(letters[1:4], 4)) Z <- im(m, xcol=1:4, yrow=1:4) # make a point pattern set.seed(42) X <- runifpoint(20, win=as.owin(Z)) # look up the image at the points of X # (a) internal ans1 <- lookup.im(Z, X$x, X$y) stopifnot(is.factor(ans1)) # (b) user level ans2 <- Z[X] stopifnot(is.factor(ans2)) # (c) turn the image into a tessellation # and apply quadratcount V <- tess(image = Z) quadratcount(X, tess=V) }) # # tests/gpclib.R # # Test whether polygon geometry works with gpclib enabled # (only testable when gpclib is installed!) # # $Revision: 1.3 $ $Date: 2011/12/05 07:29:16 $ # require(spatstat) local({ if(require(gpclib)) { spot <- spatstat.options() spatstat.options(gpclib=TRUE) example(intersect.owin) example(plot.owin) spatstat.options(spot) } }) # # tests/splitpea.R # # Check behaviour of split.ppp etc # # Thanks to Marcelino de la Cruz # # $Revision: 1.7 $ $Date: 2011/12/05 07:29:16 $ # require(spatstat) if(require(gpclib)) spatstat.options(gpclib=TRUE) local({ W <- square(8) X <- ppp(c(2.98, 4.58, 7.27, 1.61, 7.19), c(7.56, 5.29, 5.03, 0.49, 1.65), window=W) Z <- quadrats(W, 4, 4) Yall <- split(X, Z, drop=FALSE) Ydrop <- split(X, Z, drop=TRUE) P <- Yall[[1]] if(!all(inside.owin(P$x, P$y, P$window))) stop("Black hole detected when drop=FALSE") P <- Ydrop[[1]] if(!all(inside.owin(P$x, P$y, P$window))) stop("Black hole detected when drop=TRUE") Ydrop[[1]] <- P[1] split(X, Z, drop=TRUE) <- Ydrop # test NA handling Zbad <- quadrats(square(4), 2, 2) Ybdrop <- split(X, Zbad, drop=TRUE) Yball <- split(X, Zbad, drop=FALSE) # From Marcelino set.seed(1) W<- square(10) # the big window puntos<- rpoispp(0.5, win=W) data(letterR) r00 <- letterR r05 <- shift(letterR,c(0,5)) r50 <- shift(letterR,c(5,0)) r55 <- shift(letterR,c(5,5)) tessr4 <- tess(tiles=list(r00, r05,r50,r55)) puntosr4 <- split(puntos, tessr4, drop=TRUE) split(puntos, tessr4, drop=TRUE) <- puntosr4 }) # # tests/imageops.R # # $Revision: 1.5 $ $Date: 2011/12/05 07:29:16 $ # require(spatstat) local({ A <- as.im(owin()) B <- as.im(owin(c(1.1, 1.9), c(0,1))) Z <- imcov(A, B) stopifnot(abs(max(Z) - 0.8) < 0.1) }) # tests/triplets.R # test code for triplet interaction # $Revision: 1.4 $ $Date: 2012/07/12 02:43:32 $ require(spatstat) local({ fit <- ppm(redwood, ~1, Triplets(0.1)) fit suffstat(fit) # hard core (zero triangles, coefficient is NA) fit0 <- ppm(cells, ~1, Triplets(0.05)) fit0 suffstat(fit0) # bug case (1 triangle in data) fit1 <- ppm(cells, ~1, Triplets(0.15)) fit1 suffstat(fit1) }) # # tests/project.ppm.R # # $Revision: 1.3 $ $Date: 2012/10/22 03:12:08 $ # # Tests of projection mechanism # require(spatstat) local({ # a very unidentifiable model fit <- ppm(cells, ~Z, Strauss(1e-06), covariates=list(Z=0)) project.ppm(fit) # multitype fit2 <- ppm(amacrine, ~1, MultiStrauss(types=c("off", "on"), radii=matrix(1e-06, 2, 2))) project.ppm(fit2) # hybrids r0 <- min(nndist(redwood)) ra <- 1.25 * r0 rb <- 0.8 * r0 f <- ppm(redwood, ~1, Hybrid(A=Strauss(ra), B=Geyer(0.1, 2)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Strauss(rb), B=Geyer(0.1, 2)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Strauss(ra), B=Strauss(0.1)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Strauss(rb), B=Strauss(0.1)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Hardcore(rb), B=Strauss(0.1)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Hardcore(rb), B=Geyer(0.1, 2)), project=TRUE) f <- ppm(redwood, ~1, Hybrid(A=Geyer(rb, 1), B=Strauss(0.1)), project=TRUE) }) # # tests/hyperframe.R # # test "[.hyperframe" etc # # $Revision: 1.2 $ $Date: 2012/01/31 11:04:44 $ # lambda <- runif(4, min=50, max=100) X <- lapply(as.list(lambda), function(x) { rpoispp(x) }) h <- hyperframe(lambda=lambda, X=X) h$lambda2 <- lambda^2 h[, "lambda3"] <- lambda^3 h[, "Y"] <- X h[, "X"] <- lapply(X, flipxy) h[, c("X", "Y")] <- hyperframe(X=X, Y=X) # check fast code for Kest require(spatstat) local({ Kb <- Kest(cells, nlarge=0) Ku <- Kest(cells, correction="none") Kbu <- Kest(cells, correction=c("none", "border")) }) # # tests/vcovppm.R # # Check validity of vcov.ppm algorithms # # Thanks to Ege Rubak # # $Revision: 1.3 $ $Date: 2012/11/06 08:46:01 $ # require(spatstat) local({ set.seed(42) X <- rStrauss(200, .5, .05) model <- ppm(X, inter = Strauss(.05)) b <- vcov(model, generic = TRUE, algorithm = "basic") v <- vcov(model, generic = TRUE, algorithm = "vector") vc <- vcov(model, generic = TRUE, algorithm = "vectorclip") vn <- vcov(model, generic = FALSE) disagree <- function(x, y, tol=1e-7) { max(abs(x-y)) > tol } asymmetric <- function(x) { disagree(x, t(x)) } if(asymmetric(b)) stop("Non-symmetric matrix produced by vcov.ppm 'basic' algorithm") if(asymmetric(v)) stop("Non-symmetric matrix produced by vcov.ppm 'vector' algorithm") if(asymmetric(vc)) stop("Non-symmetric matrix produced by vcov.ppm 'vectorclip' algorithm") if(asymmetric(vn)) stop("Non-symmetric matrix produced by vcov.ppm Strauss algorithm") if(disagree(v, b)) stop("Disagreement between vcov.ppm algorithms 'vector' and 'basic' ") if(disagree(v, vc)) stop("Disagreement between vcov.ppm algorithms 'vector' and 'vectorclip' ") if(disagree(vn, vc)) stop("Disagreement between vcov.ppm generic and Strauss algorithms") # Geyer code xx <- c(0.7375956, 0.6851697, 0.6399788, 0.6188382) yy <- c(0.5816040, 0.6456319, 0.5150633, 0.6191592) Y <- ppp(xx, yy, window=square(1)) modelY <- ppm(Y, ~1, Geyer(0.1, 1)) b <- vcov(modelY, generic = TRUE, algorithm = "basic") v <- vcov(modelY, generic = TRUE, algorithm = "vector") vc <- vcov(modelY, generic = TRUE, algorithm = "vectorclip") if(asymmetric(b)) stop("Non-symmetric matrix produced by vcov.ppm 'basic' algorithm for Geyer model") if(asymmetric(v)) stop("Non-symmetric matrix produced by vcov.ppm 'vector' algorithm for Geyer model") if(asymmetric(vc)) stop("Non-symmetric matrix produced by vcov.ppm 'vectorclip' algorithm for Geyer model") if(disagree(v, b)) stop("Disagreement between vcov.ppm algorithms 'vector' and 'basic' for Geyer model") if(disagree(v, vc)) stop("Disagreement between vcov.ppm algorithms 'vector' and 'vectorclip' for Geyer model") }) # tests/windows.R # Tests of owin geometry code require(spatstat) local({ spatstat.options(gpclib=TRUE) # Ege Rubak spotted this problem in 1.28-1 A <- as.owin(ants) B <- dilation(A, 140) if(!is.subset.owin(A, B)) stop("is.subset.owin fails in polygonal case") spatstat.options(gpclib=FALSE) }) # # test addvar options # X <- rpoispp(function(x,y){exp(3+3*x)}) model <- ppm(X, ~y) addvar(model, "x", crosscheck=TRUE) addvar(model, "x", bw.input="quad") w <- square(0.5) addvar(model, "x", subregion=w) addvar(model, "x", subregion=w, bw.input="points") # additional test of parres X <- rpoispp(function(x,y){exp(3+x+2*x^2)}) model <- ppm(X, ~x+y) # options in parres parres(model, "x") parres(model, "x", bw.input="quad") w <- square(0.5) parres(model, "x", subregion=w) parres(model, "x", subregion=w, bw.input="quad") # check whether 'update.ppm' has messed up internals mod2 <- update(model, ~x) parres(mod2, "x") # # tests/lppstuff.R # # Tests for lpp code # # $Revision: 1.2 $ $Date: 2013/01/23 08:20:41 $ require(spatstat) local({ # check 'normalise' option in linearKinhom X <- rpoislpp(5, simplenet) fit <- lppm(X, ~x) K <- linearKinhom(X, lambda=fit, normalise=FALSE) plot(K) g <- linearpcfinhom(X, lambda=fit, normalise=FALSE) plot(g) K <- linearKinhom(X, lambda=fit, normalise=TRUE) plot(K) g <- linearpcfinhom(X, lambda=fit, normalise=TRUE) plot(g) # check empty patterns OK X <- runiflpp(0, simplenet) print(X) }) # # tests/density.R # # Test behaviour of density methods and inhomogeneous summary functions # # $Revision: 1.1 $ $Date: 2013/02/26 09:13:52 $ # require(spatstat) local({ lam <- density(redwood) K <- Kinhom(redwood, lam) lamX <- density(redwood, at="points") KX <- Kinhom(redwood, lamX) }) # # tests/slrm.R # # $Revision: 1.1 $ $Date: 2013/04/19 10:14:52 $ # # Test slrm fitting and prediction when there are NA's # require(spatstat) local({ X <- copper$SouthPoints W <- owin(poly=list(x=c(0,35,35,1),y=c(1,1,150,150))) Y <- X[W] fit <- slrm(Y ~ x+y) pred <- predict(fit) }) # tests/linalgeb.R # checks validity of linear algebra code # $Revision: 1.2 $ $Date: 2013/04/18 09:14:37 $ require(spatstat) local({ p <- 3 n <- 4 x <- array(as.numeric(1:(p * n * n)), dim=c(p, n, n)) w <- matrix(1:(n*n), n, n) y <- matrix(numeric(p * p), p, p) for(i in 1:n) for(j in (1:n)[-i]) y <- y + w[i,j] * outer(x[,i,j], x[,j,i]) z <- sumsymouter(x, w) if(!identical(y,z)) stop("sumsymouter gives incorrect result") })