tests/testsL.R

In spatstat.linnet: Linear Networks Functionality of the 'spatstat' Family

#'
#'   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)
FULLTEST <- (nchar(Sys.getenv("SPATSTAT_TEST", unset="")) > 0)
ALWAYS   <- TRUE
cat(paste("--------- Executing",
          if(FULLTEST) "** ALL **" else "**RESTRICTED** subset of",
          "test code -----------\n"))
#
# 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")
  }
})

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)
  }
  
})