R Under development (unstable) (2024-01-29 r85841 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.explore
> #'   Obtain environment variable controlling tests.
> #'
> #'   $Revision: 1.5 $ $Date: 2020/04/30 05:31:37 $
> 
> require(spatstat.explore)
Loading required package: spatstat.explore
Loading required package: spatstat.data
Loading required package: spatstat.geom
spatstat.geom 3.2-8
Loading required package: spatstat.random
spatstat.random 3.2-2
Loading required package: nlme
spatstat.explore 3.2-6
> 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/kernelstuff.R
> #'
> #'  $Revision: 1.2 $  $Date: 2023/11/05 01:49:45 $
> 
> local({
+   if(FULLTEST) {
+     #' test all cases in kernels.R
+     kernames <- c("gaussian", "rectangular", "triangular",
+                   "epanechnikov", "biweight", "cosine", "optcosine")
+     X <- rnorm(20)
+     U <- runif(20)
+     for(ker in kernames) {
+       dX <- dkernel(X, ker)
+       fX <- pkernel(X, ker)
+       qU <- qkernel(U, ker)
+       m0 <- kernel.moment(0, 0, ker)
+       m1 <- kernel.moment(1, 0, ker)
+       m2 <- kernel.moment(2, 0, ker)
+       m3 <- kernel.moment(3, 0, ker)
+     }
+   }
+ })
> 
> #'
> #'   tests/Kfuns.R
> #'
> #'   Various K and L functions and pcf
> #'
> #'   $Revision: 1.43 $  $Date: 2022/06/17 01:47:08 $
> #'
> #'   Assumes 'EveryStart.R' was run
> 
> if(FULLTEST) {
+   Cells <- cells
+   Amacrine <- amacrine
+   Redwood <- redwood
+ } else {
+   ## reduce numbers of data + dummy points
+   spatstat.options(npixel=32, ndummy.min=16)
+   Cells <- cells[c(FALSE, TRUE)]
+   Amacrine <- amacrine[c(FALSE, TRUE)]
+   Redwood <- redwood[c(FALSE, TRUE)]
+ }
> 
> 
> myfun <- function(x,y){(x+1) * y } # must be outside
> 
> local({
+   if(FULLTEST) {
+     #' supporting code
+     rmax.rule("Kscaled", owin(), 42)
+     implemented.for.K(c("border", "bord.modif", "translate", "good", "best"),
+                       "polygonal", TRUE)
+     implemented.for.K(c("border", "bord.modif", "translate", "good", "best"),
+                       "mask", TRUE)
+     implemented.for.K(c("border", "isotropic"), "mask", TRUE)
+     implemented.for.K(c("border", "isotropic"), "mask", FALSE)
+     #' shortcuts
+     D <- density(Cells)
+     K <- Kborder.engine(Cells, rmax=0.4, weights=D, ratio=TRUE)
+     K <- Knone.engine(Cells, rmax=0.4, weights=D, ratio=TRUE)
+     allcor <- c("none", "border", "bord.modif","isotropic", "translate")
+     K <- Krect.engine(Cells, rmax=0.4, ratio=TRUE, correction=allcor)
+     K <- Krect.engine(Cells, rmax=0.4, ratio=TRUE, correction=allcor,
+                       weights=D)
+     K <- Krect.engine(Cells, rmax=0.4, ratio=TRUE, correction=allcor,
+                       use.integers=FALSE)
+     #' Kest special code blocks
+     K <- Kest(Cells, var.approx=TRUE, ratio=FALSE)
+     Z <- distmap(Cells) + 1
+     Kb <- Kest(Cells, correction=c("border","bord.modif"),
+                weights=Z, ratio=TRUE)
+     Kn <- Kest(Cells, correction="none",
+                weights=Z, ratio=TRUE)
+     Knb <- Kest(Cells, correction=c("border","bord.modif","none"),
+                 weights=Z, ratio=TRUE)
+   }
+   if(ALWAYS) {
+     bigint <- 50000 # This is only "big" on a 32-bit system where
+                     # sqrt(.Machine$integer.max) = 46340.9
+     X <- runifpoint(bigint)
+     Z <- as.im(1/bigint, owin())
+     Kb <- Kest(X, correction=c("border","bord.modif"),
+                rmax=0.02, weights=Z, ratio=TRUE)
+   }
+   if(FULLTEST) {
+     Kn <- Kest(X, correction="none",
+                rmax=0.02, weights=Z, ratio=TRUE)
+     Knb <- Kest(X, correction=c("border","bord.modif","none"),
+                 rmax=0.02, weights=Z, ratio=TRUE)
+     #' pcf.ppp special code blocks
+     pr  <- pcf(Cells, ratio=TRUE, var.approx=TRUE)
+     pc  <- pcf(Cells, domain=square(0.5))
+     pcr <- pcf(Cells, domain=square(0.5), ratio=TRUE)
+     pw <- pcf(Redwood, correction="none")
+     pwr <- pcf(Redwood, correction="none", ratio=TRUE)
+     pv <- pcf(Redwood, kernel="rectangular")
+     p1 <- pcf(Redwood[1])
+     #' pcf.fv
+     K <- Kest(Redwood)
+     g <- pcf(K, method="a")
+     g <- pcf(K, method="c")
+     g <- pcf(K, method="d")
+     #' Kinhom code blocks
+     X <- rpoispp(function(x,y) { 100 * x }, 100, square(1))
+     lambda <- 100 * X$x
+     Kin <- Kinhom(X, lambda, correction=c("none", "border"))
+     lambda2 <- outer(lambda, lambda, "*")
+     Ki2 <- Kinhom(X, lambda2=lambda2, diagonal=FALSE,
+                   correction=c("translate", "isotropic"))
+   }
+   if(ALWAYS) {
+     #' edge corrections
+     rr <- rep(0.1, npoints(Cells))
+     eC <- edge.Ripley(Cells, rr)
+     eI <- edge.Ripley(Cells, rr, method="interpreted")
+     if(max(abs(eC-eI)) > 0.1)
+       stop("Ripley edge correction results do not match")
+   }
+   if(FULLTEST) {
+     a <- rmax.Ripley(square(1))
+     a <- rmax.Rigid(square(1))
+     a <- rmax.Ripley(as.polygonal(square(1)))
+     a <- rmax.Rigid(as.polygonal(square(1)))
+     a <- rmax.Ripley(letterR)
+     a <- rmax.Rigid(letterR)
+   }
+   if(ALWAYS) {
+     #' run slow code for edge correction and compare results
+     op <- spatstat.options(npixel=128)
+     X <- Redwood[c(TRUE, FALSE, FALSE, FALSE)]
+     Window(X) <- as.polygonal(Window(X))
+     Eapprox <- edge.Trans(X)
+     Eexact <- edge.Trans(X, exact=TRUE)
+     maxrelerr <- max(abs(1 - range(Eapprox/Eexact)))
+     if(maxrelerr > 0.1)
+       stop(paste("Exact and approximate algorithms for edge.Trans disagree by",
+                  paste0(round(100*maxrelerr), "%")),
+            call.=FALSE)
+     spatstat.options(op)
+   }
+ })
> 
> local({
+   if(FULLTEST) {
+     #' ----  multitype ------
+     K <- Kcross(Amacrine, correction=c("none", "bord.modif"))
+     K <- Kcross(Amacrine, correction=c("none", "bord", "bord.modif"),
+                           ratio=TRUE)
+     #' inhomogeneous multitype
+     K2 <- Kcross.inhom(Amacrine, lambdaX=densityfun(Amacrine))
+     K3 <- Kcross.inhom(Amacrine, lambdaX=density(Amacrine, at="points"))
+     K5 <- Kcross.inhom(Amacrine, correction="bord.modif")
+     #' markconnect, markcorr
+     M <- markconnect(Amacrine, "on", "off", normalise=TRUE)
+     M <- markcorr(longleaf, normalise=TRUE,
+                   correction=c("isotropic", "translate", "border", "none"))
+     M <- markcorr(longleaf, normalise=TRUE, fargs=list())
+     #' Kmark (=markcorrint)
+     X <- runifpoint(100) %mark% runif(100)
+     km <- Kmark(X, f=atan2)
+     km <- Kmark(X, f1=sin)
+     km <- Kmark(X, f="myfun")
+     aa <- Kmark(X, normalise=FALSE, returnL=FALSE)
+     aa <- Kmark(X, normalise=FALSE, returnL=TRUE)
+     aa <- Kmark(X, normalise=TRUE,  returnL=FALSE)
+     aa <- Kmark(X, normalise=TRUE,  returnL=TRUE)
+   }
+ })
> 
> local({
+   if(FULLTEST) {
+     #'    various modified K functions
+     #'
+     #'   directional K functions
+     #'
+     a <- Ksector(swedishpines,
+                  -pi/2, pi/2, units="radians",
+                  correction=c("none", "border", "bord.modif",
+                               "Ripley", "translate"),
+                  ratio=TRUE)
+     plot(a)
+     #'
+     #'   local K functions
+     #'
+     Z <- as.im(intensity(swedishpines), W=Window(swedishpines))
+     ZX <- Z[swedishpines]
+     a <- localLinhom(swedishpines, lambda=Z)
+     a <- localLinhom(swedishpines, lambda=ZX)
+     a <- localLinhom(swedishpines, lambda=Z, correction="none")
+     a <- localLinhom(swedishpines, lambda=Z, correction="translate")
+     a <- localLcross(Amacrine)
+     a <- localLcross(Amacrine, from="off", to="off")
+     a <- localKdot(Amacrine)
+     a <- localLdot(Amacrine)
+     a <- localKcross.inhom(Amacrine)
+     a <- localLcross.inhom(Amacrine)
+     Zed <- solapply(intensity(amacrine), as.im, W=Window(amacrine))
+     Lum <- evaluateCovariateAtPoints(Zed, Amacrine)
+     moff <- (marks(Amacrine) == "off")
+     a <- localLcross.inhom(Amacrine, from="off", to="on", lambdaX=Zed)
+     a <- localLcross.inhom(Amacrine, from="off", to="on", lambdaX=Lum)
+     a <- localLcross.inhom(Amacrine, from="off", to="on",
+                            lambdaFrom=Lum[moff], lambdaTo=Lum[!moff])
+     a <- localLcross.inhom(Amacrine, from="off", to="on", lambdaX=Zed,
+                            correction="none")
+     a <- localLcross.inhom(Amacrine, from="off", to="on", lambdaX=Zed,
+                            correction="translate")
+     #'
+     #' cases of resolve.lambdacross
+     #'
+     h <- resolve.lambdacross(Amacrine, moff, !moff)
+     h <- resolve.lambdacross(Amacrine, moff, !moff, lambdaX=Zed)
+     h <- resolve.lambdacross(Amacrine, moff, !moff, lambdaX=Lum)
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaI=Zed[["off"]], lambdaJ=Zed[["on"]])
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaI=Lum[moff], lambdaJ=Lum[!moff])
+     d <- densityfun(unmark(Amacrine), sigma=0.1)
+     dm <- lapply(split(Amacrine), densityfun, sigma=0.1)
+     h <- resolve.lambdacross(Amacrine, moff, !moff, lambdaX=d)
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaI=dm[["off"]], lambdaJ=dm[["on"]])
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaX=function(x,y,m){ d(x,y) })
+     #'
+     #' multitype inhomogeneous pcf
+     #'
+     g <- pcfcross.inhom(Amacrine, 
+                         lambdaI=dm[["off"]], lambdaJ=dm[["on"]])
+ 
+     #'
+     #'   lohboot code blocks
+     #'
+     Ared <- lohboot(Redwood, fun="Kest", block=TRUE,
+                     Vcorrection=TRUE, global=FALSE, correction="none")
+     Bred <- lohboot(Redwood, block=TRUE, basicboot=TRUE, global=FALSE)
+     Cred <- lohboot(Redwood, fun=Kest, block=TRUE, global=TRUE,
+                     correction="translate")
+     Dred <- lohboot(Redwood, Lest)
+     Kred <- lohboot(Redwood, Kinhom)
+     Lred <- lohboot(Redwood, Linhom)
+     gred <- lohboot(Redwood, pcfinhom, sigma=0.1)
+     #'
+     X <- runifpoint(100, letterR)
+     AX <- lohboot(X, block=TRUE, nx=7, ny=10)
+     #'    multitype
+     b <- lohboot(Amacrine, Kcross)
+     b <- lohboot(Amacrine, Lcross)
+     b <- lohboot(Amacrine, Kdot)
+     b <- lohboot(Amacrine, Ldot)
+     b <- lohboot(Amacrine, Kcross.inhom)
+     b <- lohboot(Amacrine, Lcross.inhom)
+     
+     ##  Kscaled
+     A <- Lscaled(japanesepines, renormalise=TRUE, correction="all")
+   }
+ })
>   
> local({
+   if(ALWAYS) {
+     #' From Ege, in response to a stackoverflow question.
+     #' The following example has two points separated by r = 1 with 1/4 of the
+     #' circumference outside the 10x10 window (i.e. area 100).
+     #' Thus the value of K^(r) should jump from 0 to 
+     #' 100/(2\cdot 1)\cdot ((3/4)^{-1} + (3/4)^{-1}) = 100 \cdot 4/3 = 133.333.
+     x <- c(4.5,5.5)
+     y <- c(10,10)-sqrt(2)/2
+     W <- square(10)
+     X <- ppp(x, y, W)
+     compere <- function(a, b, where, tol=1e-6) {
+       descrip <- paste("discrepancy in isotropic edge correction", where)
+       err <- as.numeric(a) - as.numeric(b)
+       maxerr <- max(abs(err))
+       blurb <- paste(descrip, "is", paste0(signif(maxerr, 4), ","), 
+                      if(maxerr > tol) "exceeding" else "within",
+                      "tolerance of", tol)
+       message(blurb)
+       if(maxerr > tol) {
+         message(paste("Discrepancies:", paste(err, collapse=", ")))
+         stop(paste("excessive", descrip), call.=FALSE)
+       }
+       invisible(TRUE)
+     }
+     ## Testing:
+     eX <- edge.Ripley(X, c(1,1))
+     compere(eX, c(4/3,4/3), "at interior point of rectangle")
+     ## Corner case:
+     Y <- X
+     Y$x <- X$x-4.5+sqrt(2)/2
+     eY <- edge.Ripley(Y, c(1,1))
+     compere(eY, c(2,4/3), "near corner of rectangle")
+     ## Invoke polygonal code
+     Z <- rotate(Y, pi/4)
+     eZdebug <- edge.Ripley(Z, c(1,1), internal=list(debug=TRUE))
+     compere(eZdebug, c(2,4/3), "at interior point of polygon (debug on)")
+     ## test validity without debugger,in case of quirks of compiler optimisation
+     eZ <- edge.Ripley(Z, c(1,1))
+     compere(eZ,      c(2,4/3), "at interior point of polygon (debug off)")
+   }
+ })
discrepancy in isotropic edge correction at interior point of rectangle is 2.22e-16, within tolerance of 1e-06
discrepancy in isotropic edge correction near corner of rectangle is 8.882e-16, within tolerance of 1e-06
/// Debug level 3 ///
------- centre[0] = (-6.071068, 7.071068) ------
        boundary distance 0.707107 
radius[0] = 1.000000
... Edge[0] = (0.000000,0.000000) to (7.071068,7.071068)
Left: det = -35.857864
Right: det = -171.715729
Top: det = -34143.145751
Finished cutting; ncut = 0
contrib = -0.000000
... Edge[1] = (7.071068,7.071068) to (0.000000,14.142136)
Left: det = -171.715729
Right: det = -35.857864
Top: det = -34143.145751
Finished cutting; ncut = 0
contrib = 0.000000
... Edge[2] = (0.000000,14.142136) to (-7.071068,7.071068)
Left: det = -35.857864
Right: det = 0.000000
	det = 0
Top: det = 200.000000
	det > 0
	hits + segment: t = 1.000000, theta = 3.141593
	hits - segment: t = 0.858579, theta = 1.570796
Finished cutting; ncut = 2
theta[0] = 1.570796
theta[1] = 3.141593
Interval 0, width 1.570796:inside
Interval 1, width 1.570796:outside
Interval 2, width 3.141593:inside
contrib = 4.712389
... Edge[3] = (-7.071068,7.071068) to (0.000000,0.000000)
Left: det = 0.000000
	det = 0
Right: det = -35.857864
Top: det = 200.000000
	det > 0
	hits + segment: t = 0.141421, theta = -1.570796
	hits - segment: t = 0.000000, theta = 3.141593
Finished cutting; ncut = 2
theta[0] = 3.141593
theta[1] = 4.712389
Interval 0, width 3.141593:outside
Interval 1, width 1.570796:inside
Interval 2, width 1.570796:outside
contrib = -1.570796

Total = 3.141593 = 0.500000 * (2 * pi)
------- centre[1] = (-5.363961, 7.778175) ------
        boundary distance 0.707107 
radius[0] = 1.000000
... Edge[0] = (0.000000,0.000000) to (7.071068,7.071068)
Left: det = -27.772078
Right: det = -153.629942
Top: det = -34143.145751
Finished cutting; ncut = 0
contrib = -0.000000
... Edge[1] = (7.071068,7.071068) to (0.000000,14.142136)
Left: det = -153.629942
Right: det = -27.772078
Top: det = -27108.831175
Finished cutting; ncut = 0
contrib = 0.000000
... Edge[2] = (0.000000,14.142136) to (-7.071068,7.071068)
Left: det = -27.772078
Right: det = -1.914214
Top: det = 200.000000
	det > 0
	hits + segment: t = 0.900000, theta = 3.141593
	hits - segment: t = 0.758579, theta = 1.570796
Finished cutting; ncut = 2
theta[0] = 1.570796
theta[1] = 3.141593
Interval 0, width 1.570796:inside
Interval 1, width 1.570796:outside
Interval 2, width 3.141593:inside
contrib = 4.712389
... Edge[3] = (-7.071068,7.071068) to (0.000000,0.000000)
Left: det = -1.914214
Right: det = -27.772078
Top: det = -765.685425
Finished cutting; ncut = 0
contrib = -0.000000

Total = 4.712389 = 0.750000 * (2 * pi)
discrepancy in isotropic edge correction at interior point of polygon (debug on) is 7.55e-15, within tolerance of 1e-06
discrepancy in isotropic edge correction at interior point of polygon (debug off) is 7.55e-15, within tolerance of 1e-06
> 
> 
> 
> reset.spatstat.options()
> 
> proc.time()
   user  system elapsed 
   2.01    0.29    2.28