R Under development (unstable) (2024-09-18 r87177 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #' > #' Header for all (concatenated) test files > #' > #' Require spatstat.linnet > #' Obtain environment variable controlling tests. > #' > #' $Revision: 1.5 $ $Date: 2020/04/30 05:31:37 $ > > require(spatstat.linnet) Loading required package: spatstat.linnet Loading required package: spatstat.data Loading required package: spatstat.univar spatstat.univar 3.0-1 Loading required package: spatstat.geom spatstat.geom 3.3-3 Loading required package: spatstat.random spatstat.random 3.3-2 Loading required package: spatstat.explore Loading required package: nlme spatstat.explore 3.3-2 Loading required package: spatstat.model Loading required package: rpart spatstat.model 3.3-2 spatstat.linnet 3.2-2 > FULLTEST <- (nchar(Sys.getenv("SPATSTAT_TEST", unset="")) > 0) > ALWAYS <- TRUE > cat(paste("--------- Executing", + if(FULLTEST) "** ALL **" else "**RESTRICTED** subset of", + "test code -----------\n")) --------- Executing **RESTRICTED** subset of test code ----------- > # > # tests/lppstuff.R > # > # Tests for lpp code > # > # $Revision: 1.72 $ $Date: 2024/06/05 09:11:58 $ > > > local({ + if(ALWAYS) { + #' make test data + Xsimple <- runiflpp(6, simplenet) + Xchic <- chicago[c(TRUE, FALSE)] + Xspid <- spiders[c(TRUE, FALSE)] + Xdend <- dendrite[seq_len(npoints(dendrite)) %% 8 == 0] + } + + if(FULLTEST) { + #' lpp class support + Xone <- Xsimple %mark% runif(6) + Xtwo <- Xsimple %mark% data.frame(a=1:6, b=runif(6)) + print(summary(Xone)) + print(summary(Xtwo)) + plot(Xsimple, show.window=TRUE) + plot(Xone) + plot(Xtwo, do.several=FALSE) + #' geometry etc + rotate(Xsimple, pi/3, centre=c(0.2,0.3)) + superimpose.lpp(L=simplenet) + W <- Window(Xsimple) + #' cut.lpp + tes <- lineardirichlet(Xsimple[1:4]) + f <- as.linfun(tes) + Z <- as.linim(f) + cut(Xsimple, tes) + cut(Xsimple, f) + cut(Xsimple, Z) + + #' check 'normalise' option in linearKinhom + fit <- lppm(Xsimple ~x) + K <- linearKinhom(Xsimple, lambda=fit, normalise=FALSE, leaveoneout=FALSE) + plot(K) + g <- linearpcfinhom(Xsimple, lambda=fit, normalise=FALSE, leaveoneout=FALSE) + plot(g) + K <- linearKinhom(Xsimple, lambda=fit, normalise=TRUE, leaveoneout=FALSE) + plot(K) + g <- linearpcfinhom(Xsimple, lambda=fit, normalise=TRUE, leaveoneout=FALSE) + plot(g) + ## other code blocks + K <- linearKinhom(Xsimple, lambda=fit, correction="none", ratio=TRUE) + g <- linearpcf(Xsimple, correction="none", ratio=TRUE) + g1 <- linearpcf(Xsimple[1], ratio=TRUE) + K1 <- linearKcross(dendrite[1], "thin", "thin", ratio=TRUE) + ## check empty patterns OK + X0 <- runiflpp(0, simplenet) + print(X0) + g <- linearpcf(X0, ratio=TRUE) + } + + ## nearest neighbour distances + 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) } + + if(ALWAYS) { + nn <- nndist(Xspid) + nnP <- f(pairdist(Xspid), 1) + if(any(abs(nn - nnP) > eps)) + stop("nndist.lpp does not agree with pairdist.lpp") + + nw <- nnwhich(Xspid) + nwP <- g(pairdist(Xspid), 1) + if(any(nw != nwP)) + stop("nnwhich.lpp does not agree with pairdist") + } + + if(FULLTEST) { + #' code blocks in nndist.lpp/nnwhich.lpp + #' non-sparse network, interpreted code + Ad <- nndist(Xspid, method="interpreted") + Aw <- nnwhich(Xspid, method="interpreted") + #' sparse network, older C code + opa <- spatstat.options(Cnndistlpp=FALSE) + Bd <- nndist(dendrite) + Bw <- nnwhich(dendrite) + spatstat.options(opa) + #' undefined nearest neighbours + Ed <- nndist(Xspid[1:3], k=1:3) + Ew <- nnwhich(Xspid[1:3], k=1:3) + + #' trivial cases in nncross.lpp + a <- nncross(runiflpp(0, simplenet), runiflpp(1, simplenet), + what="which", format="list")$which + a <- nncross(runiflpp(0, simplenet), runiflpp(1, simplenet), + what="dist", format="list")$dist + } + + if(ALWAYS) { + #' compare algorithms + ZZ <- split(Xchic) + XX <- ZZ$damage + YY <- ZZ$assault + op <- spatstat.options(Cnncrosslpp=FALSE) + a <- nncross(XX, YY) + spatstat.options(Cnncrosslpp=TRUE) + b <- nncross(XX, YY) + if(any(a$which != b$which)) + stop("Inconsistent values of nncross.lpp()$which from different C code") + if(max(abs(a$dist - b$dist)) > eps) + stop("Inconsistent values of nncross.lpp()$dist from different C code") + + spatstat.options(Cnncrosslpp=TRUE) + b2 <- nncross(XX, YY, k=1:2, what="which") + if(any(b2$which.1 != b$which)) + stop("inconsistent values of nncross.lpp()$which from k=1:2 and k=1") + a2 <- nncross(XX, YY, k=1:2, what="dist") + if(max(abs(a2$dist.1 - a$dist)) > eps) + stop("Inconsistent values of nncross.lpp()$dist from k=1:2 and k=1") + + spatstat.options(Cnncrosslpp=TRUE) + ii <- seq_len(npoints(XX)) + w1 <- nnwhich(XX) + w2 <- nncross(XX, XX, iX=ii, iY=ii, what="which") + w3 <- nncross(XX, XX, iX=ii, iY=ii, what="which", method="interpreted") + if(any(w1 != w2)) + stop("nnwhich.lpp disagrees with nncross.lpp(iX, iY)") + if(any(w2 != w3)) + stop("Different results for nncross.lpp(iX, iY, 'which') using R and C") + d1 <- nndist(XX) + d2 <- nncross(XX, XX, iX=ii, iY=ii, what="dist") + d3 <- nncross(XX, XX, iX=ii, iY=ii, what="dist", method="interpreted") + if(max(abs(d1-d2)) > eps) + stop("nndist.lpp disagrees with nncross.lpp(iX, iY)") + if(max(abs(d2-d3)) > eps) + stop("Different results for nncross.lpp(iX, iY, 'dist') using R and C") + spatstat.options(Cnncrosslpp=FALSE) + w4 <- nncross(XX, XX, iX=ii, iY=ii, what="which") + d4 <- nncross(XX, XX, iX=ii, iY=ii, what="dist") + if(any(w2 != w4)) + stop("Different results for nncross.lpp(iX, iY, 'which') fast and slow C") + if(max(abs(d2-d4)) > eps) + stop("Different results for nncross.lpp(iX, iY, 'dist') fast and slow C") + spatstat.options(Cnncrosslpp=TRUE) + } + + reset.spatstat.options() + + if(FULLTEST) { + ## test handling marginal cases + xyd <- nncross(XX, YY[1]) + A <- runiflpp(5, simplenet) + B <- runiflpp(2, simplenet) + aaa <- nncross(A,B,k=3:5) #' all undefined + aaa <- nncross(A,B,k=1:4) #' some undefined + spatstat.options(Cnncrosslpp=FALSE) + aaa <- nncross(A,B,k=3:5) + aaa <- nncross(A,B,k=1:4) + bbb <- nncross(B,A, iX=1:2, iY=1:5) # another code block + spatstat.options(Cnncrosslpp=TRUE) + reset.spatstat.options() + } + + if(FULLTEST) { + ## as.linnet.psp (Suman's example) + Lines <- as.data.frame(as.psp(simplenet)) + newseg <- c(Lines[1,1:2], Lines[10,3:4]) + Lines <- rbind(Lines, newseg) + Y <- as.psp(Lines, window=Window(simplenet)) + marks(Y) <- c(3, 4, 5, 5, 3, 4, 5, 5,5, 5,1) + Z <- as.linnet(Y) # can crash if marks don't match segments + + ## Test linnet surgery code + RL <- joinVertices(simplenet, matrix(c(2,3), ncol=2)) # redundant edge + RZ <- joinVertices(Z, matrix(c(2,3), ncol=2), marks=6) # redundant edge + JZ <- joinVertices(Z, matrix(c(2,7), ncol=2), marks=6) # new edge + set.seed(42) + X <- runiflpp(30, simplenet) + V <- runiflpp(30, simplenet) + XV <- insertVertices(X, V) + validate.lpp.coords(XV, context="calculated by insertVertices") + X0 <- insertVertices(X, x=numeric(0), y=numeric(0)) + ## vertices on boundary of new window + LL <- simplenet[boundingbox(vertices(simplenet))] + + ## thinNetwork for lpp and linim + Xthin <- thinNetwork(X, retainedges=1:3) + D <- density(X, 0.1) + Dthin <- thinNetwork(D, retainedges=1:3) + + ## Test [.lpp internal data + B <- owin(c(0.1,0.7),c(0.19,0.6)) + XB <- X[B] + validate.lpp.coords(XB, context="returned by [.lpp") + } + + ## Tests related to linearK, etc + testcountends <- function(X, r=100, s=1) { + if(s != 1) { + X <- rescale(X, s) + r <- r/s + } + L <- as.linnet(X) + n1 <- countends(L, X[1], r) + n2 <- npoints(lineardisc(L, X[1], r, plotit=FALSE)$endpoints) + if(n1 != n2) + stop(paste("Incorrect result from countends:", + n1, "!=", n2, + spatstat.utils::paren(paste("scale=", 1/s))), + call.=FALSE) + } + + if(ALWAYS) { + ## original scale + XU <- unmark(Xchic) + testcountends(XU) + ## finer scale + testcountends(XU, s=1000) + #' disconnected + L <- thinNetwork(simplenet, retainedges = -c(3,8)) + S <- as.psp(L) + x <- midpoints.psp(S)[1] + len <- lengths_psp(S)[1] + A <- lineardisc(L, x, len, plotit=FALSE) #involves many segments of network + B <- lineardisc(L, x, len/5, plotit=FALSE) # involves one segment of network + op <- spatstat.options(Ccountends=FALSE) + A <- lineardisc(L, x, len, plotit=FALSE) + B <- lineardisc(L, x, len/5, plotit=FALSE) + reset.spatstat.options() + + ## Test algorithms for boundingradius.linnet + L <- as.linnet(chicago, sparse=TRUE) + L$boundingradius <- NULL # artificially remove + opa <- spatstat.options(Clinearradius=FALSE) + bR <- as.linnet(L, sparse=FALSE)$boundingradius + spatstat.options(Clinearradius=TRUE) + bC <- as.linnet(L, sparse=FALSE)$boundingradius + spatstat.options(opa) + if(abs(bR-bC) > 0.001 * (bR+bC)/2) + stop("Disagreement between R and C algorithms for boundingradius.linnet", + call.=FALSE) + } + + if(FULLTEST) { + ## linnet things + is.connected(as.linnet(dendrite)) + zik <- rescale(Xchic, 39.37/12, "m") + Simon <- simplenet + unitname(Simon) <- list("metre", "metres", 0.5) + b <- rescale(Simon) + ds <- density(simplenet, 0.05) + } + + if(ALWAYS) { + ## invoke dist2dpath + LS <- as.linnet(simplenet, sparse=TRUE) + LF <- as.linnet(LS, sparse=FALSE) + ## direct call dist2dpath + d <- simplenet$dpath + d[!simplenet$m] <- Inf + diag(d) <- 0 + dd <- dist2dpath(d, method="interpreted") + ra <- range(dd - simplenet$dpath) + if(max(abs(ra)) > sqrt(.Machine$double.eps)) + stop("dist2dpath gives different answers in C and R code") + } + + if(FULLTEST) { + ## integral.linim with missing entries + xcoord <- linfun(function(x,y,seg,tp) { x }, domain(Xchic)) + xcoord <- as.linim(xcoord, dimyx=32) + integral(xcoord) + + ## Math.linim + range(xcoord) + stopifnot(is.linim(sqrt(xcoord))) + stopifnot(is.linim(-xcoord)) + stopifnot(is.linim(xcoord + xcoord)) + stopifnot(is.linim(xcoord/3)) + + ## options to plot.linim + plot(xcoord, legend=FALSE) + plot(xcoord, leg.side="top") + plot(xcoord, style="width", leg.side="bottom") + + ## as.linim.linim + xxcc <- as.linim(xcoord) + xxcceps <- as.linim(xcoord, eps=15) + xxccdel <- as.linim(xcoord, delta=30) + df1 <- attr(xxcc, "df") + df2 <- attr(xxccdel, "df") + df3 <- resampleNetworkDataFrame(df1, df2) + + ## linim with complex values + Zc <- as.im(function(x,y){(x-y) + x * 1i}, Frame(simplenet)) + Fc <- linim(simplenet, Zc) + print(Fc) + print(summary(Fc)) + + ## linim with df provided + Z <- as.im(function(x,y) {x-y}, Frame(simplenet)) + X <- linim(simplenet, Z) + df <- attr(X, "df") + XX <- linim(simplenet, Z, df=df) + dfwithout <- df[, colnames(df) != "values"] + XXX <- linim(simplenet, Z, df=dfwithout) + plot(XXX, zlim=c(-1,1)) + plot(XXX, legend=FALSE) + plot(XXX, leg.side="bottom") + + ## lpp with multiple columns of marks + M <- Xchic + marks(M) <- cbind(type=marks(M), data.frame(distnearest=nndist(M))) + plot(M, main="") + summary(M) + MM <- cut(M) + + #' other cases + CC <- cut(Xchic) + nd <- nndist(Xspid) + SX <- cut(Xspid %mark% nd, breaks=3) + SX <- cut(Xspid, nd, breaks=c(0,100,200,Inf), include.lowest=TRUE) + + ## linequad + Y <- Xsimple %mark% factor(rep(c("A", "B"), 3)) + aX <- linequad(Xsimple) + aY <- linequad(Y) + aXR <- linequad(Xsimple, random=TRUE) + aYR <- linequad(Y, random=TRUE) + P <- as.ppp(Xsimple) + S <- as.psp(domain(Xsimple)) + d <- linequad(P, S) + oop <- spatstat.options(Clinequad=FALSE) + bX <- linequad(Xsimple) + spatstat.options(oop) + + ## other internal utilities + df <- pointsAlongNetwork(simplenet, 0.05) + X <- as.ppp(df[,c("x", "y")], W=Frame(simplenet)) + A <- local2lpp(simplenet, seg=df$seg, tp=df$tp, X=X, df.only=FALSE) + + ## mark-mark scatterplot uses pairdist + X <- runiflpp(20, simplenet) %mark% runif(20) + markmarkscatter(X, 0.2) + markmarkscatter(X[FALSE], 0.1) + + ## tree branches + ## make a simple tree + m <- simplenet$m + m[8,10] <- m[10,8] <- FALSE + L <- linnet(vertices(simplenet), m) + tb <- treebranchlabels(L, 1) + X <- runiflpp(30, L) + ## delete branch B + XminusB <- deletebranch(X, "b", tb) + ## extract branch B + XB <- extractbranch(X, "b", tb) + + ## cases of lintess() + A <- lintess(simplenet) # argument df missing + S <- as.psp(simplenet) + ns <- nsegments(S) + df <- data.frame(seg=1:ns, t0=0, t1=1, tile=letters[1:ns]) + M <- data.frame(len=lengths_psp(S), ang=angles.psp(S)) + V <- lintess(simplenet, df, marks=M) + + ## methods for class lintess + U <- unmark(V) + U <- unstack(V) + print(summary(V)) + W <- Window(V) + plot(V, style="image") + plot(V, style="width") + + ## linear tessellations infrastructure + nX <- 100 + nY <- 20 + X <- runiflpp(nX, simplenet) + Y <- runiflpp(nY, simplenet) + tes <- divide.linnet(Y) + cX <- coords(X) + iI <- lineartileindex(cX$seg, cX$tp, tes, method="interpreted") + iC <- lineartileindex(cX$seg, cX$tp, tes, method="C") + iE <- lineartileindex(cX$seg, cX$tp, tes, method="encode") + if(!identical(iI,iC)) + stop("conflicting results from lineartileindex (interpreted vs C)") + if(!identical(iI,iE)) + stop("conflicting results from lineartileindex (interpreted vs encoded)") + iA <- as.linfun(tes)(X) + if(!identical(iI, iA)) + stop("disagreement between as.linfun.lintess and lineartileindex") + + ## intersection of lintess + X <- divide.linnet(runiflpp(4, simplenet)) + Y <- divide.linnet(runiflpp(3, simplenet)) + nX <- summary(X)$nt + nY <- summary(Y)$nt + marks(X) <- factor(letters[1L + (seq_len(nX) %% 2)]) + marks(Y) <- runif(nY) + Zmm <- intersect.lintess(X,Y) + Zum <- intersect.lintess(unmark(X),Y) + Zmu <- intersect.lintess(X,unmark(Y)) + } + + if(FULLTEST) { + #' handling by 'solist', 'unstack', 'plot.solist' etc + L <- simplenet + X <- runiflpp(5, L) %mark% cbind(a=1:5, b=letters[1:5]) + ns <- nsegments(L) + df <- data.frame(seg=1:ns, t0=0, t1=1, tile=letters[1:ns]) + S <- lintess(L, df) + f <- as.linfun(S) + g <- as.linfun(S, values=seq_len(nsegments(L))) + V <- as.linim(f) + Z <- as.linim(g) + shebang <- solist(L=L, X=X, S=S, f=f, g=g, V=V, Z=Z) + plot(shebang) + plot(shebang, valuesAreColours=FALSE) + kapow <- unstack(shebang) + plot(kapow) + } + + #' kernel density estimation + X <- Xsimple[1:5] + w <- runif(npoints(X)) + if(ALWAYS) { + #' density.lpp -> densityHeat, includes platform-dependent code + D <- density(X, 0.05) + D <- density(X, 0.05, weights=w) + } + if(FULLTEST) { + #' code blocks + D <- density(X[FALSE], 0.05) + D <- density(X, Inf) + D <- density(X, 0.05, finespacing=TRUE) # } + D <- density(X, 0.05, eps=0.008) # } code blocks in resolve.heat.steps + D <- density(X, 0.05, dimyx=256) # } + } + if(ALWAYS) { + DX <- density(X, 0.05, at="points", fastmethod="a", debug=TRUE) + DX <- density(X, 0.05, at="points", fast=FALSE, debug=TRUE) + ## densityfun.lpp, code blocks + ff <- densityfun(X, 0.05, nsigma=Inf) + } + if(FULLTEST) { + #' disconnected network + L <- thinNetwork(simplenet, retainedges=-c(3,5)) + Y <- runiflpp(5, L) + D <- density(Y, 0.05) + D <- density(Y, 0.05, weights=w) + D <- density(Y, Inf) + g <- flatdensityfunlpp(Y, disconnect=FALSE) + a <- flatdensityatpointslpp(Y, disconnect=FALSE) + #' density.lpp -> densityEqualSplit + D <- density(X, 0.05, kernel="e", weights=w) + D <- density(X, Inf, kernel="e") + D <- density(Y, 0.05, kernel="e", weights=w) + D <- density(Y, Inf, kernel="e") + } + if(ALWAYS) { + #' density.lpp -> densityQuick.lpp + D <- density(X, 0.05, distance="e", weights=runif(5)) + } + if(FULLTEST) { + D <- density(X, Inf, distance="e") + D <- density(Y, 0.05, distance="e", weights=runif(5)) + D <- density(Y, Inf, distance="e") + ## example from Andrea Gilardi + ## with duplicated points, and points at a vertex + A <- runiflpp(5, simplenet) + V <- lpp(vertices(simplenet)[2:4], simplenet) + B <- superimpose(A, A[3], V) + D <- density(B, sigma=0.2, distance="e", at="points") + if(length(as.numeric(D)) != npoints(B)) + stop("density.lpp(at='points') returns wrong number of values") + G <- superimpose(U=runiflpp(4, simplenet), B=B) + b <- bw.relrisklpp(G, distance="e", nh=2) + } + if(FULLTEST) { + #' densityVoronoi.lpp and related code + densityVoronoi(X, f=0) + densityVoronoi(X, f=1e-8) + densityVoronoi(X, f=1) + densityVoronoi(X[FALSE], f=0.5) + XX <- X[rep(1:npoints(X), 4)] + densityVoronoi(XX, f=0.99999, nrep=5) + densityVoronoi(Y, f=0) + densityVoronoi(Y, f=1e-8) + densityVoronoi(Y, f=1) + densityVoronoi(Y[FALSE], f=0.5) + #' bandwidth selection + bw.voronoi(X, nrep=4, prob=c(0.2, 0.4, 0.6)) + bw.lppl(X) + sug <- seq(0.1, 0.3, by=0.05) + bw.lppl(X, sigma=sug) + bw.lppl(X, sigma=sug, shortcut=FALSE) + bw.lppl(X, sigma=sug, distance="path", shortcut=FALSE) + } + if(ALWAYS) { + #' inhomogeneous K and g + SD <- split(Xdend) + DT <- density(SD[["thin"]], 100, distance="euclidean") + DS <- density(SD[["stubby"]], 100, distance="euclidean") + Kii <- linearKcross.inhom(dendrite, "thin", "thin", DT, DT) + Kij <- linearKcross.inhom(dendrite, "thin", "stubby", DT, DS, + correction="none", ratio=TRUE) + gii <- linearpcfcross.inhom(dendrite, "thin", "thin", DT, DT) + gij <- linearpcfcross.inhom(dendrite, "thin", "stubby", DT, DS, + correction="none", ratio=TRUE) + } + if(FULLTEST) { + Kx <- linearKcross.inhom(dendrite[1], "thin", "stubby", DT, DS) + gx <- linearpcfcross.inhom(dendrite[1], "thin", "stubby", DT, DS) + } + + if(FULLTEST) { + ## infrastructure for density.lpp + L <- domain(Xchic) + A <- resolve.heat.steps(100, L=L, dx=1) + B <- resolve.heat.steps(100, L=L, dt=0.2) + C <- resolve.heat.steps(100, L=L, niter=1e5) + C <- resolve.heat.steps(100, L=L, niter=1e5, nsave=3) + C <- resolve.heat.steps(100, L=L, niter=1e5, nsave=Inf) + D <- resolve.heat.steps(100, L=L, dx=1, dt=0.2) + E <- resolve.heat.steps(500, L=L, dt=0.5, iterMax=2e5) + A <- resolve.heat.steps(1, L=simplenet, dt=2, dx=0.05) + B <- resolve.heat.steps(1, L=simplenet) + C <- resolve.heat.steps(1, L=simplenet, finespacing=FALSE) + D <- resolve.heat.steps(1, seglengths=rep(0.7, 5), maxdegree=4, AMbound=7) + ## vertices on boundary + L <- simplenet + L <- L[boundingbox(vertices(L))] + X <- as.lpp(0.41259, 0.6024, L=L) + D <- density(X, 0.1) + } + + if(FULLTEST) { + #' multitype + #' density.splitppx + Z <- split(Xchic)[1:3] + D <- density(Z, 7) + #' relrisk.lpp for 3 types + U <- do.call(superimpose, Z) + A <- relrisk(U, 20, control=1, relative=TRUE) + A <- relrisk(U, 20, control=1, relative=TRUE, at="points") + #' relrisk.lpp for 2 types + V <- do.call(superimpose, Z[1:2]) + A <- relrisk(V, 20, control=2, casecontrol=FALSE) # issues warning + B <- relrisk(V, 20, control="assault", case=2, relative=TRUE) + D <- relrisk(V, 20, control=1, case="burglary", relative=TRUE, at="points") + #' bw.relrisklpp + b1 <- bw.relrisklpp(V, hmax=3, fudge=0.01, + method="leastsquares", reference="sigma") + b2 <- bw.relrisklpp(V, hmax=3, fudge=0.01, + method="KelsallDiggle", reference="uniform") + } + + if(FULLTEST) { + #' pairs.linim + Z <- density(Xsimple, 0.5, distance="euclidean") + pairs(solist(Z)) + pairs(solist(A=Z)) + U <- density(as.ppp(Xsimple), 0.5) + pairs(solist(U, Z)) + pairs(solist(Z, U)) + } + + if(FULLTEST) { + #' complex-valued functions and images + f <- function(x,y,seg,tp) { x + y * 1i } + g <- linfun(f, simplenet) + h <- as.linim(g) + plot(Re(h)) + plot(h) + plot(g) + integral(h) + integral(g) + } + + ## 'lixellate' - involves C code + if(ALWAYS) { + ## Cases where no subdivision occurs + P <- runiflpp(4, simplenet) + A <- lixellate(P, nsplit=1) + B <- lixellate(P, eps=2) + ## bug in 'lixellate' (Jakob Gulddahl Rasmussen) + X <- ppp(c(0,1), c(0,0), owin()) + L <- linnet(X, edges = matrix(1:2, ncol=2)) + Y <- lpp(X, L) + ## The left end point is OK + lixellate(Y[1], nsplit=30) + d <- density(Y[1], .1) + ## The right end point gave an error + lixellate(Y[2], nsplit=30) + d <- density(Y[2], .1) + } + + if(FULLTEST) { + ## make some bad data and repair it + L <- simplenet + ## reverse edges + a <- L$from[c(FALSE,TRUE)] + L$from[c(FALSE,TRUE)] <- L$to[c(FALSE,TRUE)] + L$to[c(FALSE,TRUE)] <- a + ## duplicate edges + ns <- nsegments(L) + ii <- c(seq_len(ns), 2) + L$from <- L$from[ii] + L$to <- L$to[ii] + L$lines <- L$lines[ii] + ## + Y <- repairNetwork(L) + ## add points + X <- runiflpp(4, L) + Z <- repairNetwork(X) + } + + if(FULLTEST) { + #' random generation bugs and code blocks + A <- runiflpp(5, simplenet, nsim=2) + D <- density(A[[1]], 0.3) + B <- rpoislpp(D, nsim=2) + stopifnot(is.multitype(rlpp(c(10,5), list(a=D,b=D)))) + stopifnot(is.multitype(rlpp(5, list(a=D,b=D)))) + stopifnot(is.multitype(rlpp(c(10,5), D))) + } + + ## rhohat.lppm + if(FULLTEST) { + fut <- lppm(Xspid ~ 1) + rx <- rhohat(fut, "x") + Z <- linfun(function(x,y,seg,tp) { x }, domain(Xspid)) + rZ <- rhohat(fut, Z) + U <- predict(rx) + U <- predict(rZ) + Y <- simulate(rx) + Y <- simulate(rZ) + futm <- lppm(Xchic ~ x + marks) + ry <- rhohat(futm, "y") + U <- predict(ry) + Y <- simulate(ry) + } + + if(ALWAYS) { + #' new code for pairdist.lpp + X <- runiflpp(10, simplenet) + dC <- pairdist(X, method="C") + dI <- pairdist(X, method="interpreted") + if(max(abs(dC - dI)) > 0.01) + stop("pairdist.lpp: disagreement between C and R") + XS <- as.lpp(X, sparse=TRUE) + dS <- pairdist(XS) + if(max(abs(dC - dS)) > 0.01) + stop("pairdist.lpp: disagreement between sparse and non-sparse C") + dU <- pairdist(XS, method="testsymm") # secret option + if(max(abs(dS - dU)) > 0.01) + stop("pairdist.lpp: disagreement between symmetric and asymmetric") + + #' new code for crossdist.lpp + #' non-sparse + v <- split(Xchic) + X <- v$cartheft + Y <- v$burglary + dC <- crossdist(X, Y, method="C") + dI <- crossdist(X, Y, method="interpreted") + if(max(abs(dC - dI)) > 0.01) + stop("crossdist.lpp (non-sparse): disagreement between C and R") + #' convert to sparse + Chuck <- as.lpp(Xchic, sparse=TRUE) + V <- split(Chuck) + X <- V$cartheft + Y <- V$burglary + dS <- crossdist(X, Y) + if(max(abs(dS - dC)) > 0.01) + stop("crossdist.lpp: disagreement between sparse and non-sparse") + #' disconnected + L <- thinNetwork(simplenet, retainedges = -c(3,8)) + X <- runiflpp(20, L) + Y <- runiflpp(15, L) + d <- crossdist(X,Y) + + ## example where the path is very long (covers a previous bug) + X <- dendrite[c(349,563)] + a <- pairdist(X, method="testsymm") + if(max(abs(a-t(a))) > 0.01) + stop("pairdist.lpp: asymmetry for long distances") + b <- pairdist(as.lpp(X, sparse=FALSE)) + if(max(abs(a-b)) > 0.01) + stop("pairdist.lpp: disagreement sparse vs non-sparse for long distance") + } + }) sigma = 0.05 Network: total length = 2.90485162043096 number of edges = 10 average nonzero edge length = 0.290485162043096 shortest nonzero edge length = 0.150146561419501 Determine dx by default rule Mean Nonzero Edge Length/ 30 = 0.00968283873476987 Total Network Length/ 100 = 0.0290485162043096 Min Nonzero Edge Length/3 = 0.0500488538065003 Minimum of the above = 0.00968283873476987 Fine spacing rule: dx = 0.00968283873476987 Rounded total number of lixels = 300 Adjusted dx = 0.00968283873476987 Applying full set of constraints dt (determined by all constraints) = 1.27242139521974e-05 Number of iterations (determined from dt) = 98 Updated dt = 1.27551020408163e-05 alpha = 0.136043732774807 Final values: Time step dt = 1.27551020408163e-05 Sample spacing dx = 0.00968283873476987 Number of iterations niter = 98 Number of states saved nsave = 1 Subdividing network ...Done. New network: Linear network with 299 vertices and 299 lines Vertices are marked Enclosing window: rectangle = [0, 1] x [0, 1] units Constructing update matrix A ..Done. alpha = 0.136043732774807 Building initial state ..Done. Running iterative solver ..Done. Mapping results to spatial location ..Done. sigma = 0.05 Network: total length = 2.90485162043096 number of edges = 10 average nonzero edge length = 0.290485162043096 shortest nonzero edge length = 0.150146561419501 Determine dx by default rule Mean Nonzero Edge Length/ 30 = 0.00968283873476987 Total Network Length/ 100 = 0.0290485162043096 Min Nonzero Edge Length/3 = 0.0500488538065003 Minimum of the above = 0.00968283873476987 Fine spacing rule: dx = 0.00968283873476987 Rounded total number of lixels = 300 Adjusted dx = 0.00968283873476987 Applying full set of constraints dt (determined by all constraints) = 1.27242139521974e-05 Number of iterations (determined from dt) = 98 Updated dt = 1.27551020408163e-05 alpha = 0.136043732774807 Final values: Time step dt = 1.27551020408163e-05 Sample spacing dx = 0.00968283873476987 Number of iterations niter = 98 Number of states saved nsave = 1 Subdividing network ...Done. New network: Linear network with 299 vertices and 299 lines Vertices are marked Enclosing window: rectangle = [0, 1] x [0, 1] units Constructing update matrix A ..Done. alpha = 0.136043732774807 Building initial state ..Done. finenet: Linear network with 299 vertices and 299 lines Vertices are marked Enclosing window: rectangle = [0, 1] x [0, 1] units lixelmap: int [1:10] 261 262 283 284 231 232 163 164 218 219 lixelweight: num [1:10] 1.02 102.25 14.22 89.05 83.3 ... Processing 5 points ... 1, 2, 3, 4, 5. > > reset.spatstat.options() > > #' > #' lppmodels.R > #' > #' Tests of lppm and class support > #' > #' $Revision: 1.2 $ $Date: 2020/06/12 00:24:09 $ > #' > > > local({ + if(ALWAYS) { + fit0 <- lppm(spiders) + fit1 <- lppm(spiders ~ x) + summary(fit0) + summary(fit1) + pseudoR2(fit0) + pseudoR2(fit1) + } + if(FULLTEST) { + fit2 <- lppm(chicago ~ x+y) + summary(fit2) + pseudoR2(fit2) + } + if(ALWAYS) { + X <- runiflpp(10, simplenet) + Z <- distfun(runiflpp(10, simplenet)) + fit3 <- lppm(X ~ Z) + summary(fit3) + pseudoR2(fit3) + } + + Window(fit1) + + if(ALWAYS) a <- model.images(fit0) + if(FULLTEST) { + a <- model.images(fit1) + a <- model.images(fit2) + a <- model.images(fit3) + } + if(ALWAYS) + b <- model.matrix(fit0) + if(FULLTEST) { + b <- model.matrix(fit1) + b <- model.matrix(fit2) + b <- model.matrix(fit3) + } + + if(ALWAYS) + is.multitype(fit0) + if(FULLTEST) { + is.multitype(fit1) + is.multitype(fit2) + is.multitype(fit3) + } + + if(ALWAYS) fit0e <- emend(fit0) + if(FULLTEST) { + fit1e <- emend(fit1) + fit2e <- emend(fit2) + fit3e <- emend(fit3) + } + + #' fundamental utilities: + #' spatialCovariateEvidence (was: evalCovar) + ycoord <- function(x,y) { y } + if(ALWAYS) YS <- as.linim(ycoord, L=domain(spiders)) + if(FULLTEST) YC <- as.linim(ycoord, L=domain(chicago)) + + if(ALWAYS) aT <- spatialCovariateEvidence(fit1, YS, interpolate=TRUE) + if(FULLTEST) { + aF <- spatialCovariateEvidence(fit1, YS, interpolate=FALSE) + dT <- spatialCovariateEvidence(fit1, ycoord, interpolate=TRUE) + dF <- spatialCovariateEvidence(fit1, ycoord, interpolate=FALSE) + bT <- spatialCovariateEvidence(fit2, YC, interpolate=TRUE) + bF <- spatialCovariateEvidence(fit2, YC, interpolate=FALSE) + cT <- spatialCovariateEvidence(fit2, ycoord, interpolate=TRUE) + cF <- spatialCovariateEvidence(fit2, ycoord, interpolate=FALSE) + } + + }) Point process model on linear network Fitted to point pattern dataset 'spiders' Internal fit: Point process model Fitted to data: linequad(spiders, eps = NULL, nd = 1000, random = FALSE) Fitting method: maximum likelihood (Berman-Turner approximation) Model was fitted using glm() Algorithm converged Call: ppm.quad(Q = linequad(spiders, eps = NULL, nd = 1000, random = FALSE), method = "mpl", forcefit = TRUE) Edge correction: "border" [border correction distance r = 0 ] -------------------------------------------------------------------------------- Quadrature scheme (Berman-Turner) = data + dummy + weights Data pattern: Planar point pattern: 48 points Average intensity 3.79e-05 points per square mm Window: rectangle = [0, 1125] x [0, 1125] mm Window area = 1265620 square mm Unit of length: 1 mm Dummy quadrature points: Equally spaced along each segment at spacing eps = 20.22 mm Original parameter nd = 1000 Planar point pattern: 1198 points Average intensity 0.000947 points per square mm Window: rectangle = [0, 1125] x [0, 1125] mm Window area = 1265620 square mm Unit of length: 1 mm Quadrature weights: Counting weights based on segment length All weights: range: [3.12, 20.2] total: 20200 Weights on data points: range: [4.78, 10.1] total: 466 Weights on dummy points: range: [3.12, 20.2] total: 19800 -------------------------------------------------------------------------------- FITTED : Stationary Poisson process ---- Intensity: ---- Uniform intensity: [1] 0.002374034 Estimate S.E. CI95.lo CI95.hi Ztest Zval (Intercept) -6.043165 0.1443376 -6.326061 -5.760268 *** -41.86827 ----------- gory details ----- Fitted regular parameters (theta): (Intercept) -6.043165 Fitted exp(theta): (Intercept) 0.002374034 Original data: spiders Domain: Linear network with 156 vertices and 203 lines Total length 20218.75 mm Maximum vertex degree: 3 [Non-sparse matrix representation] Network is connected Diameter: 2575 mm Bounding radius: 1306.25 mm Numerical tolerance: 0.00625 mm Enclosing window: rectangle = [0, 1125] x [0, 1125] mm Point process model on linear network Fitted to point pattern dataset 'spiders' Internal fit: Point process model Fitted to data: linequad(spiders, eps = NULL, nd = 1000, random = FALSE) Fitting method: maximum likelihood (Berman-Turner approximation) Model was fitted using glm() Algorithm converged Call: ppm.quad(Q = linequad(spiders, eps = NULL, nd = 1000, random = FALSE), trend = ~x, interaction = NULL, data = NULL, method = "mpl", forcefit = TRUE) Edge correction: "border" [border correction distance r = 0 ] -------------------------------------------------------------------------------- Quadrature scheme (Berman-Turner) = data + dummy + weights Data pattern: Planar point pattern: 48 points Average intensity 3.79e-05 points per square mm Window: rectangle = [0, 1125] x [0, 1125] mm Window area = 1265620 square mm Unit of length: 1 mm Dummy quadrature points: Equally spaced along each segment at spacing eps = 20.22 mm Original parameter nd = 1000 Planar point pattern: 1198 points Average intensity 0.000947 points per square mm Window: rectangle = [0, 1125] x [0, 1125] mm Window area = 1265620 square mm Unit of length: 1 mm Quadrature weights: Counting weights based on segment length All weights: range: [3.12, 20.2] total: 20200 Weights on data points: range: [4.78, 10.1] total: 466 Weights on dummy points: range: [3.12, 20.2] total: 19800 -------------------------------------------------------------------------------- FITTED : Nonstationary Poisson process ---- Intensity: ---- Log intensity: ~x Fitted trend coefficients: (Intercept) x -6.4845961867 0.0007423366 Estimate S.E. CI95.lo CI95.hi Ztest (Intercept) -6.4845961867 0.3179208001 -7.1077095048 -5.861482869 *** x 0.0007423366 0.0004467129 -0.0001332047 0.001617878 Zval (Intercept) -20.396892 x 1.661775 ----------- gory details ----- Fitted regular parameters (theta): (Intercept) x -6.4845961867 0.0007423366 Fitted exp(theta): (Intercept) x 0.001526777 1.000742612 Original data: spiders Domain: Linear network with 156 vertices and 203 lines Total length 20218.75 mm Maximum vertex degree: 3 [Non-sparse matrix representation] Network is connected Diameter: 2575 mm Bounding radius: 1306.25 mm Numerical tolerance: 0.00625 mm Enclosing window: rectangle = [0, 1125] x [0, 1125] mm Point process model on linear network Fitted to point pattern dataset 'X' Internal fit: Point process model Fitted to data: linequad(X, eps = NULL, nd = 1000, random = FALSE) Fitting method: maximum likelihood (Berman-Turner approximation) Model was fitted using glm() Algorithm converged Call: ppm.quad(Q = linequad(X, eps = NULL, nd = 1000, random = FALSE), trend = ~Z, interaction = NULL, data = NULL, method = "mpl", forcefit = TRUE) Edge correction: "border" [border correction distance r = 0 ] -------------------------------------------------------------------------------- Quadrature scheme (Berman-Turner) = data + dummy + weights Data pattern: Planar point pattern: 10 points Average intensity 10 points per square unit Window: rectangle = [0, 1] x [0, 1] units Window area = 1 square unit Dummy quadrature points: Equally spaced along each segment at spacing eps = 0.002905 units Original parameter nd = 1000 Planar point pattern: 1010 points Average intensity 1010 points per square unit Window: rectangle = [0, 1] x [0, 1] units Window area = 1 square unit Quadrature weights: Counting weights based on segment length All weights: range: [0.000763, 0.0029] total: 2.9 Weights on data points: range: [0.00145, 0.00145] total: 0.0145 Weights on dummy points: range: [0.000763, 0.0029] total: 2.89 -------------------------------------------------------------------------------- FITTED : Nonstationary Poisson process ---- Intensity: ---- Log intensity: ~Z Model depends on external covariate 'Z' Covariates provided: Z: function(...) Fitted trend coefficients: (Intercept) Z 1.519107 -2.980157 Estimate S.E. CI95.lo CI95.hi Ztest Zval (Intercept) 1.519107 0.5126746 0.5142837 2.523931 ** 2.9631025 Z -2.980157 4.6104558 -12.0164841 6.056170 -0.6463909 ----------- gory details ----- Fitted regular parameters (theta): (Intercept) Z 1.519107 -2.980157 Fitted exp(theta): (Intercept) Z 4.56814625 0.05078487 Original data: X Domain: Linear network with 10 vertices and 10 lines Total length 2.904852 units Maximum vertex degree: 4 [Non-sparse matrix representation] Network is connected Diameter: 1.248943 units Bounding radius: 0.7793609 units Numerical tolerance: 0.0001501466 units Enclosing window: rectangle = [0, 1] x [0, 1] units > > > proc.time() user system elapsed 6.68 1.46 8.10