R Under development (unstable) (2024-11-18 r87347 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.model
> #'   Obtain environment variable controlling tests.
> #'
> #'   $Revision: 1.5 $ $Date: 2020/04/30 05:31:37 $
> 
> require(spatstat.model)
Loading required package: spatstat.model
Loading required package: spatstat.data
Loading required package: spatstat.univar
spatstat.univar 3.1-1
Loading required package: spatstat.geom
spatstat.geom 3.3-4
Loading required package: spatstat.random
spatstat.random 3.3-2
Loading required package: spatstat.explore
Loading required package: nlme
spatstat.explore 3.3-3
Loading required package: rpart
spatstat.model 3.3-3
> 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/kppm.R
> #
> # $Revision: 1.40 $ $Date: 2024/02/26 05:44:02 $
> #
> # Test functionality of kppm that once depended on RandomFields
> # Test update.kppm for old style kppm objects
> 
> if(!FULLTEST)
+   spatstat.options(npixel=32, ndummy.min=16)
> 
> local({
+ 
+  fit <- kppm(redwood ~1, "Thomas") # sic
+  fitx <- kppm(redwood ~x, "Thomas", verbose=TRUE)
+  if(FULLTEST) {
+    fitx <- update(fit, ~ . + x)
+    fitM <- update(fit, clusters="MatClust")
+    fitC <- update(fit, cells)
+    fitCx <- update(fit, cells ~ x)
+    #'
+    Wsub <- owin(c(0, 0.5), c(-0.5, 0))
+    Zsub <- (bdist.pixels(Window(redwood)) > 0.1)
+    fitWsub <- kppm(redwood ~1, "Thomas", subset=Wsub)
+    fitZsub <- kppm(redwood ~1, "Thomas", subset=Zsub)
+    fitWsub
+  
+    #' various methods
+    ff <- as.fv(fitx)
+    uu <- unitname(fitx)
+    unitname(fitCx) <- "furlong"
+    mo <- model.images(fitCx)
+    p <- psib(fit)
+    px <- psib(fitx)
+  }
+  if(ALWAYS) {
+    Y <- simulate(fitx, seed=42, saveLambda=TRUE)[[1]]
+  }
+ 
+  if(FULLTEST) {
+    #' vcov.kppm different algorithms
+    vc  <- vcov(fitx)
+    vc2 <- vcov(fitx, fast=TRUE)
+    vc3 <- vcov(fitx, fast=TRUE, splitup=TRUE)
+    vc4 <- vcov(fitx,            splitup=TRUE)
+ 
+    ## other code blocks
+    a <- varcount(fitx, function(x,y){x+1}) # always positive
+    a <- varcount(fitx, function(x,y){y-1}) # always negative
+    a <- varcount(fitx, function(x,y){x+y}) # positive or negative
+ 
+    #' improve.kppm
+    fitI <- update(fit, improve.type="quasi")
+    fitxI <- update(fitx, improve.type="quasi")
+    fitxIs <- update(fitx, improve.type="quasi", fast=FALSE) 
+    #' vcov.kppm
+    vcI <- vcov(fitxI)
+  }
+ 
+   ## plot.kppm including predict.kppm
+  if(ALWAYS) {
+    fitMC <- kppm(redwood ~ x, "Thomas")
+    plot(fitMC)
+  }
+  if(FULLTEST) {
+    fitCL <- kppm(redwood ~ x, "Thomas", method="c")
+    fitPA <- kppm(redwood ~ x, "Thomas", method="p")
+    plot(fitCL)
+    plot(fitPA)
+ 
+    ## fit with composite likelihood method [thanks to Abdollah Jalilian]
+    fut <- kppm(redwood ~ x, "VarGamma", method="clik2", nu.ker=-3/8)
+    kfut <- as.fv(fut)
+  }
+  
+  if(ALWAYS) {
+    fit0 <- kppm(redwood ~1, "LGCP")
+    is.poisson(fit0)
+    Y0 <- simulate(fit0, saveLambda=TRUE)[[1]]
+    stopifnot(is.ppp(Y0))
+    p0 <- psib(fit0) # issues a warning
+ 
+    if(FULLTEST) {
+      ## fit LGCP using K function: slow
+      fit1 <- kppm(redwood ~x, "LGCP",
+                   covmodel=list(model="matern", nu=0.3),
+                   control=list(maxit=3))
+      Y1 <- simulate(fit1, saveLambda=TRUE)[[1]]
+      stopifnot(is.ppp(Y1))
+    }
+ 
+    ## fit LGCP using pcf
+    fit1p <- kppm(redwood ~x, "LGCP",
+                  covmodel=list(model="matern", nu=0.3),
+                  statistic="pcf")
+    Y1p <- simulate(fit1p, saveLambda=TRUE)[[1]]
+    stopifnot(is.ppp(Y1p))
+ 
+    ## .. and using different fitting methods
+    if(FULLTEST) {
+      fit1pClik <- update(fit1p, method="clik")
+      fit1pPalm <- update(fit1p, method="palm")
+    }
+ 
+    ## shortcut evaluation of pcf
+    ##  (the code being tested is in spatstat.random::clusterinfo.R)
+    if(FULLTEST) {
+      putSpatstatVariable("RFshortcut", TRUE)
+      fitGshort <- kppm(redwood ~ 1, "LGCP", covmodel=list(model="gauss"))
+      fitSshort <- kppm(redwood ~ 1, "LGCP", covmodel=list(model="stable", alpha=1))
+      putSpatstatVariable("RFshortcut", FALSE)
+      fitGlong <-  kppm(redwood ~ 1, "LGCP", covmodel=list(model="gauss"))
+      fitSlong <-  kppm(redwood ~ 1, "LGCP", covmodel=list(model="stable", alpha=1))
+      discrepG <- unlist(parameters(fitGshort)) - unlist(parameters(fitGlong))
+      discrepS <- unlist(parameters(fitSshort)) - unlist(parameters(fitSlong))
+      print(discrepG)
+      print(discrepS)
+      if(max(abs(discrepG) > 0.01))
+        stop("Discrepancy in short-cut fitting of Gaussian LGCP")
+      if(max(abs(discrepS) > 0.01))
+        stop("Discrepancy in short-cut fitting of stable LGCP")
+    }
+    
+    ## image covariate (a different code block) 
+    xx <- as.im(function(x,y) x, Window(redwood))
+    fit1xx <- update(fit1p, . ~ xx, data=solist(xx=xx))
+    Y1xx <- simulate(fit1xx, saveLambda=TRUE)[[1]]
+    stopifnot(is.ppp(Y1xx))
+    if(FULLTEST) {
+      fit1xxVG <- update(fit1xx, clusters="VarGamma", nu=-1/4)
+      Y1xxVG <- simulate(fit1xxVG, saveLambda=TRUE)[[1]]
+      stopifnot(is.ppp(Y1xxVG))
+    }
+    fit1xxLG <- update(fit1xx, clusters="LGCP",
+                       covmodel=list(model="matern", nu=0.3),
+                       statistic="pcf")
+    Y1xxLG <- simulate(fit1xxLG, saveLambda=TRUE, drop=TRUE)
+    stopifnot(is.ppp(Y1xxLG))
+    
+    # ... and Abdollah's code
+    if(FULLTEST) {
+      fit2 <- kppm(redwood ~x, cluster="Cauchy", statistic="K")
+      Y2 <- simulate(fit2, saveLambda=TRUE)[[1]]
+      stopifnot(is.ppp(Y2))
+    }
+  }
+ 
+ })
Fitting cluster model
Retrieved cluster model information
Algorithm parameters:
$rmax
NULL

$q
[1] 0.25

$p
[1] 2

$rmin
NULL

Using point pattern data
Starting parameters:
       kappa       sigma2 
62.000000000  0.006173033 
Calculating summary function...Done.
Starting minimum contrast fit
Returned from minimum contrast fit
Returning from clusterfit
Warning messages:
1: Internal error: fvlabels truncated the function name 
2: Internal error: fvlabels truncated the function name 
3: Internal error: fvlabels truncated the function name 
4: In psib.kppm(fit0) : The model is not a cluster process
5: The value of the empirical function 'pcf' for r= 0 was Inf. Range of r values was reset to [0.00048828125, 0.25] 
6: The value of the empirical function 'pcf' for r= 0 was Inf. Range of r values was reset to [0.00048828125, 0.25] 
7: The value of the empirical function 'pcf' for r= 0 was Inf. Range of r values was reset to [0.00048828125, 0.25] 
> 
> if(FULLTEST) {
+ local({
+   #' various code blocks
+   fut <- kppm(redwood, ~x)
+   fet <- update(fut, redwood3)
+   fot <- update(fut, trend=~y)
+   fit <- kppm(redwoodfull ~ x)
+   Y <- simulate(fit, window=redwoodfull.extra$regionII, saveLambda=TRUE)
+   gut <- improve.kppm(fit, type="wclik1")
+   gut <- improve.kppm(fit, vcov=TRUE, fast.vcov=TRUE, save.internals=TRUE)
+   hut <- kppm(redwood ~ x, method="clik", weightfun=NULL)
+   hut <- kppm(redwood ~ x, method="palm", weightfun=NULL)
+   mut <- kppm(redwood)
+   nut <- update(mut, Y)
+ })
+ }
> 
> if(FULLTEST) {
+ local({
+   #' minimum contrast code
+   K <- Kest(redwood)
+   a <- matclust.estK(K)
+   a <- thomas.estK(K)
+   a <- cauchy.estK(K)
+   a <- vargamma.estK(K)
+   a <- lgcp.estK(K)
+ 
+   print(a)
+   u <- unitname(a)
+   
+   g <- pcf(redwood)
+   a <- matclust.estpcf(g)
+   a <- thomas.estpcf(g)
+   a <- cauchy.estpcf(g)
+   a <- vargamma.estpcf(g)
+   a <- lgcp.estpcf(g)
+ 
+   #' auxiliary functions
+   b <- resolve.vargamma.shape(nu.pcf=1.5)
+   Z <- clusterfield("Thomas", kappa=1, scale=0.2)
+   
+   aa <- NULL
+   aa <- accumulateStatus(simpleMessage("Woof"), aa)
+   aa <- accumulateStatus(simpleMessage("Sit"),  aa)
+   aa <- accumulateStatus(simpleMessage("Woof"), aa)
+   printStatusList(aa)
+ 
+   RMIN <- 0.01
+   fit <- kppm(redwood ~ 1, ctrl=list(rmin=RMIN,q=1/2))
+   if(fit$Fit$mcfit$ctrl$rmin != RMIN)
+     stop("kppm did not handle parameter 'rmin' in argument 'ctrl' ")
+   fit <- kppm(redwood ~ 1, ctrl=list(rmin=0,q=1/2), rmin=RMIN)
+   if(fit$Fit$mcfit$ctrl$rmin != RMIN)
+     stop("kppm did not handle parameter 'rmin' in argument 'ctrl'")
+ 
+   RMIN <- 2
+   fit <- dppm(swedishpines~1, dppGauss(), ctrl=list(rmin=RMIN,q=1))
+   if(fit$Fit$mcfit$ctrl$rmin != RMIN)
+     stop("dppm did not handle parameter 'rmin' in argument 'ctrl'")
+   fit <- dppm(swedishpines~1, dppGauss(), ctrl=list(rmin=0,q=1), rmin=RMIN)
+   if(fit$Fit$mcfit$ctrl$rmin != RMIN)
+     stop("dppm did not handle argument 'rmin'")
+ })
+ }
> 
> 
> 
> if(FULLTEST) {
+ local({
+   #' cover a few code blocks
+   fut <- kppm(redwood ~ x, method="clik")
+   print(summary(fut))
+   a <- residuals(fut)
+   fut2 <- kppm(redwood ~ x, "LGCP", method="palm")
+   print(summary(fut2))
+   b <- residuals(fut2)
+   #'
+   po <- ppm(redwood ~ 1)
+   A <- kppmComLik(redwood, Xname="redwood", po=po, clusters="Thomas",
+                   statistic="pcf", statargs=list(), control=list(),
+                   weightfun=NULL, rmax=0.1)
+   A <- kppmPalmLik(redwood, Xname="redwood", po=po, clusters="Thomas",
+                    statistic="pcf", statargs=list(), control=list(),
+                    weightfun=NULL, rmax=0.1)
+ })
+ }
> 
> reset.spatstat.options()
>   
> #'
> #'   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)]
+ }
> 
> 
> if(FULLTEST) {
+   local({
+     #' code blocks using fitted model to determine intensity
+     #' Kinhom
+     X <- rpoispp(function(x,y) { 100 * x }, 100, square(1))
+     fut <- ppm(X ~ x)
+     Kio <- Kinhom(X, fut, update=FALSE)
+     Kiu <- Kinhom(X, fut, update=TRUE, diagonal=FALSE)
+     fit <- ppm(Amacrine ~ marks)
+     #' lohboot Linhom
+     Zred <- predict(ppm(Redwood ~ x+y))
+     Lred <- lohboot(Redwood, Linhom, lambda=Zred)
+     #' Kmulti.inhom
+     K1 <- Kcross.inhom(Amacrine, lambdaX=fit)
+     On <- split(Amacrine)$on
+     Off <- split(Amacrine)$off
+     K4 <- Kcross.inhom(Amacrine, lambdaI=ppm(On), lambdaJ=ppm(Off))
+     #' local K functions
+     fut <- ppm(swedishpines ~ polynom(x,y,2))
+     Z <- predict(fut)
+     Lam <- fitted(fut, dataonly=TRUE)
+     a <- localLinhom(swedishpines, lambda=fut)
+     a <- localLinhom(swedishpines, lambda=Z)
+     a <- localLinhom(swedishpines, lambda=Lam)
+     a <- localLinhom(swedishpines, lambda=Z, correction="none")
+     a <- localLinhom(swedishpines, lambda=Z, correction="translate")
+     #' local cross K functions
+     fat <- ppm(Amacrine ~ x * marks)
+     Zed <- predict(fat)
+     Lum <- fitted(fat, dataonly=TRUE)
+     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", lambdaX=fat)
+     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, lambdaX=fat)
+     h <- resolve.lambdacross(Amacrine, moff, !moff, lambdaX=fat, update=FALSE)
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaI=fat, lambdaJ=fat)
+     h <- resolve.lambdacross(Amacrine, moff, !moff,
+                               lambdaI=fat, lambdaJ=fat,
+                               update=FALSE)
+     #' lohboot
+     b <- lohboot(Amacrine, Lcross.inhom, from="off", to="on", lambdaX=Zed)
+     b <- lohboot(Amacrine, Lcross.inhom, from="off", to="on", lambdaX=Lum)
+     b <- lohboot(Amacrine, Lcross.inhom, from="off", to="on", lambdaX=fat)
+     b <- lohboot(Amacrine, Lcross.inhom, from="off", to="on", 
+                  lambdaFrom=Lum[moff], lambdaTo=Lum[!moff])
+     #'
+     #'  residual K functions etc
+     #'
+     rco <- compareFit(Cells, Kcom,
+                       interaction=anylist(P=Poisson(), S=Strauss(0.08)),
+                       same="trans", different="tcom")
+     fit <- ppm(Cells ~ x, Strauss(0.07))
+     K <- Kcom(Cells, model=fit, restrict=TRUE)
+   })
+ }
> 
> reset.spatstat.options()
> 
> proc.time()
   user  system elapsed 
   3.76    0.45    4.18