R Under development (unstable) (2024-07-08 r86886 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.0-0 Loading required package: spatstat.geom spatstat.geom 3.3-2 Loading required package: spatstat.random spatstat.random 3.3-1 Loading required package: spatstat.explore Loading required package: nlme spatstat.explore 3.3-1 Loading required package: rpart spatstat.model 3.3-1 > 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 4.15 0.31 4.46