Skip to content

Commit

Permalink
cleaned up the code
Browse files Browse the repository at this point in the history
  • Loading branch information
EricMarcon committed Nov 7, 2024
1 parent 2d6ac4b commit 8f6c2ad
Show file tree
Hide file tree
Showing 6 changed files with 172 additions and 64 deletions.
1 change: 1 addition & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
## Improvements

- avoided accessing ppp objects' marks directly.
- cleaned up the code


# dbmss 2.9-2
Expand Down
34 changes: 25 additions & 9 deletions R/DEnvelope.R
Original file line number Diff line number Diff line change
@@ -1,22 +1,38 @@
DEnvelope <-
function(X, r = NULL, NumberOfSimulations = 100, Alpha = 0.05,
Cases, Controls, Intertype = FALSE, Global = FALSE, verbose = interactive()) {
DEnvelope <- function(
X,
r = NULL,
NumberOfSimulations = 100,
Alpha = 0.05,
Cases,
Controls,
Intertype = FALSE,
Global = FALSE,
verbose = interactive()) {

CheckdbmssArguments()

# The only null hypothesis is random labelling (equivalently, random location)
SimulatedPP <- expression(rRandomLocation(X, CheckArguments = FALSE))

# local envelope, keep extreme values for lo and hi (nrank=1)
Envelope <- envelope(X, fun=Dhat, nsim=NumberOfSimulations, nrank=1,
r=r, Cases=Cases, Controls=Controls, Intertype=Intertype,
CheckArguments = FALSE,
simulate=SimulatedPP, verbose=verbose, savefuns=TRUE
)
Envelope <- envelope(
X,
fun = Dhat,
nsim = NumberOfSimulations,
nrank = 1,
r = r,
Cases = Cases,
Controls = Controls,
Intertype = Intertype,
CheckArguments = FALSE,
simulate = SimulatedPP,
verbose = verbose,
savefuns = TRUE
)
attr(Envelope, "einfo")$H0 <- "Random Location"

# Calculate confidence intervals
Envelope <- FillEnvelope(Envelope, Alpha, Global)
# Return the envelope
return (Envelope)
return(Envelope)
}
61 changes: 49 additions & 12 deletions R/Dhat.R
Original file line number Diff line number Diff line change
@@ -1,33 +1,70 @@
Dhat <-
function(X, r = NULL, Cases, Controls = NULL, Intertype = FALSE, CheckArguments = TRUE) {
Dhat <- function(
X,
r = NULL,
Cases,
Controls = NULL,
Intertype = FALSE,
CheckArguments = TRUE) {

if (CheckArguments) {
CheckdbmssArguments()
}
# K of cases.
KCases <- Khat(X, r, Cases, Cases, CheckArguments = FALSE)
KCases <- Khat(
X,
r = r,
ReferenceType = Cases,
NeighborType = Cases,
CheckArguments = FALSE
)
# Default controls are all points except cases. Reserved name is "CoNtRoLs_"
Y <- X
if (is.null(Controls)) {
Controls <- "CoNtRoLs_"
if (Controls %in% levels(marks(Y)$PointType)) stop("A point type is named 'CoNtRoLs_'. It must be changed to use the 'Controls = NULL' option of Dhat.")
if (Controls %in% levels(marks(Y)$PointType)) {
stop("A point type is named 'CoNtRoLs_'. It must be changed to use the 'Controls = NULL' option of Dhat.")
}
levels(marks(Y)$PointType) <- c(levels(marks(Y)$PointType), Controls)
marks(Y)$PointType[marks(Y)$PointType != Cases] <- Controls
}
# K of controls. r must be those of cases.
if (Intertype) {
KControls <- Khat(Y, KCases$r, Cases, Controls, CheckArguments = FALSE)
KControls <- Khat(
Y,
r = KCases$r,
ReferenceType = Cases,
NeighborType = Controls,
CheckArguments = FALSE
)
} else {
KControls <- Khat(Y, KCases$r, Controls, Controls, CheckArguments = FALSE)
KControls <- Khat(
Y,
r = KCases$r,
ReferenceType = Controls,
NeighborType = Controls,
CheckArguments = FALSE
)
}
# Calculate the difference (a difference between fv's yields a dataframe)
Dvalues <- KCases-KControls
DEstimate <- cbind(as.data.frame(KCases)[1], as.data.frame(Dvalues)[2:3])
Dvalues <- KCases - KControls
DEstimate <- cbind(
as.data.frame(KCases)[1],
as.data.frame(Dvalues)[2:3]
)

# Return the values of D(r)
D <- fv(DEstimate, argu="r",
ylab=quote(D(r)), valu=attr(KCases, "valu"), fmla=attr(KCases, "fmla"),
alim=attr(KCases, "alim"), labl=c("r", "%s[theo](r)", "hat(%s)[iso](r)"),
desc=attr(KCases, "desc"), unitname=attr(KCases, "unitname"), fname="D")
D <- fv(
DEstimate,
argu="r",
ylab = quote(D(r)),
valu = attr(KCases, "valu"),
fmla = attr(KCases, "fmla"),
alim = attr(KCases, "alim"),
labl = c("r", "%s[theo](r)", "hat(%s)[iso](r)"),
desc = attr(KCases, "desc"),
unitname = attr(KCases, "unitname"),
fname = "D"
)
fvnames(D, ".") <- colnames(DEstimate)[-1]
return (D)
}
47 changes: 30 additions & 17 deletions R/Dtable.R
Original file line number Diff line number Diff line change
@@ -1,16 +1,22 @@
Dtable <-
function(Dmatrix, PointType = NULL, PointWeight = NULL)
{
Dtable <- function(
Dmatrix,
PointType = NULL,
PointWeight = NULL) {

# Dmatrix must be a distance matrix
if (is.matrix(Dmatrix)) {
if (nrow(Dmatrix) != ncol(Dmatrix))
if (nrow(Dmatrix) != ncol(Dmatrix)) {
stop("Dmatrix should be a square matrix.")
if (any(is.na(Dmatrix)))
}
if (any(is.na(Dmatrix))) {
stop("NAs are not allowed in the distance matrix.")
if (any(Dmatrix < 0))
}
if (any(Dmatrix < 0)) {
stop("negative values are not allowed in the distance matrix.")
if (any(diag(Dmatrix) > 0))
}
if (any(diag(Dmatrix) > 0)) {
stop("diagonal values of the distance matrix must be 0.")
}
} else {
stop("Dmatrix must be a matrix.")
}
Expand All @@ -32,31 +38,38 @@ function(Dmatrix, PointType = NULL, PointWeight = NULL)
}

# Check PointType
if (any(is.null(PointType)))
if (any(is.null(PointType))) {
stop("NULL values are not allowed in the point types.")
if (any(is.na(PointType)))
}
if (any(is.na(PointType))) {
stop("NAs are not allowed in the point types.")
if (length(PointType) != nrow(Dmatrix))
}
if (length(PointType) != nrow(Dmatrix)) {
stop("The vector of point types must have the same size as Dmatrix.")
}
PointType <- as.factor(PointType)

# Get PointWeight
if (is.null(PointWeight)) {
PointWeight <- rep(1, length(PointType))
} else {
if (any(is.na(PointWeight)))
if (any(is.na(PointWeight))) {
stop("NAs are not allowed in the point weights.")
if (any(PointWeight <= 0))
}
if (any(PointWeight <= 0)) {
stop("Point weights must be strictly positive.")
if (length(PointWeight) != nrow(Dmatrix))
}
if (length(PointWeight) != nrow(Dmatrix)) {
stop("The vector of point weights must have the same size as Dmatrix.")
}
}

# Build the object
Dt <- list(Dmatrix=Dmatrix,
n=nrow(Dmatrix),
marks=list(PointType=PointType, PointWeight=PointWeight)
)
Dt <- list(
Dmatrix = Dmatrix,
n = nrow(Dmatrix),
marks = list(PointType = PointType, PointWeight = PointWeight)
)
class(Dt) <- "Dtable"
return (Dt)
}
23 changes: 16 additions & 7 deletions R/FillEnvelope.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
FillEnvelope <-
function(Envelope, Alpha, Global) {
FillEnvelope <- function(
Envelope,
Alpha,
Global) {

# Sims contains simulated values (one line per simulation)
Sims <- attr(Envelope, "simfuns")
Expand All @@ -11,10 +13,17 @@ function(Envelope, Alpha, Global) {
attr(Envelope, "desc")[5] <- "upper global envelope of %s from simulations"
attr(Envelope, "einfo")$global <- TRUE
} else {
Env <- apply(as.data.frame(Sims)[, -1], 1, stats::quantile, probs=c(Alpha/2, 1-Alpha/2), na.rm=TRUE)
# quantile returns 0 when a value is NA. Force NA to appear in the result: add NA or 0.
Env <- t(t(Env) + apply(as.data.frame(Sims), 1, sum)*0)
attr(Envelope, "einfo")$nrank <- attr(Envelope, "einfo")$nsim*Alpha/2
Env <- apply(
as.data.frame(Sims)[, -1],
MARGIN = 1,
FUN = stats::quantile,
probs = c(Alpha / 2, 1 - Alpha / 2),
na.rm = TRUE
)
# quantile returns 0 when a value is NA.
# Force NA to appear in the result: add NA or 0.
Env <- t(t(Env) + apply(as.data.frame(Sims), 1, sum) * 0)
attr(Envelope, "einfo")$nrank <- attr(Envelope, "einfo")$nsim * Alpha / 2
}
# Computed values are injected into the envelope object
Envelope$lo <- Env[1,]
Expand All @@ -23,4 +32,4 @@ function(Envelope, Alpha, Global) {
attr(Envelope, "einfo")$Alpha <- Alpha
class(Envelope) <- c("dbmssEnvelope", class(Envelope))
return(Envelope)
}
}
70 changes: 51 additions & 19 deletions R/envelope.Dtable.R
Original file line number Diff line number Diff line change
@@ -1,17 +1,29 @@
envelope.Dtable <-
function (Y, fun = Kest, nsim = 99, nrank = 1, ..., funargs = list(),
funYargs = funargs, simulate = NULL, verbose = TRUE, savefuns = FALSE,
Yname = NULL, envir.simul = NULL)
{
envelope.Dtable <- function (
Y,
fun = Kest,
nsim = 99,
nrank = 1,
...,
funargs = list(),
funYargs = funargs,
simulate = NULL,
verbose = TRUE,
savefuns = FALSE,
Yname = NULL,
envir.simul = NULL) {

# Environment information sent to envelopeEngine
cl <- spatstat.utils::short.deparse(sys.call())
if (is.null(Yname))
Yname <- spatstat.utils::short.deparse(substitute(Y))
if (is.null(fun))
fun <- Kest
envir.user <- if (!is.null(envir.simul))
if (is.null(fun)) {
fun <- Kest
}
envir.user <- if (!is.null(envir.simul)) {
envir.simul
else parent.frame()
} else {
parent.frame()
}
envir.here <- sys.frame(sys.nframe())

if (is.null(simulate)) {
Expand All @@ -22,14 +34,34 @@ function (Y, fun = Kest, nsim = 99, nrank = 1, ..., funargs = list(),
}

# Run the simulations
envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim = nsim,
nrank = nrank, ..., funargs = funargs, funYargs = funYargs,
verbose = verbose, clipdata = FALSE, transform = NULL,
global = FALSE, ginterval = NULL, use.theory = NULL,
alternative = c("two.sided", "less", "greater"), scale = NULL,
clamp = FALSE,
savefuns = savefuns, savepatterns = FALSE, nsim2 = nsim,
VARIANCE = FALSE, nSD = 2, Yname = Yname, maxnerr = nsim,
cl = cl, envir.user = envir.user, do.pwrong = FALSE,
foreignclass = "Dtable")
envelopeEngine(
X = X,
fun = fun,
simul = simrecipe,
nsim = nsim,
nrank = nrank,
...,
funargs = funargs,
funYargs = funYargs,
verbose = verbose,
clipdata = FALSE,
transform = NULL,
global = FALSE,
ginterval = NULL,
use.theory = NULL,
alternative = c("two.sided", "less", "greater"),
scale = NULL,
clamp = FALSE,
savefuns = savefuns,
savepatterns = FALSE,
nsim2 = nsim,
VARIANCE = FALSE,
nSD = 2,
Yname = Yname,
maxnerr = nsim,
cl = cl,
envir.user = envir.user,
do.pwrong = FALSE,
foreignclass = "Dtable"
)
}

0 comments on commit 8f6c2ad

Please sign in to comment.