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 8, 2024
1 parent 4423110 commit d4f3504
Show file tree
Hide file tree
Showing 18 changed files with 534 additions and 288 deletions.
31 changes: 20 additions & 11 deletions R/GlobalEnvelope.R
Original file line number Diff line number Diff line change
@@ -1,34 +1,43 @@
GlobalEnvelope <-
function(Simulations, Alpha) {
GlobalEnvelope <- function(
Simulations,
Alpha) {

verifyclass(Simulations, "fv")

# Initialize
Simulations <- as.data.frame(Simulations)[, -1] # Eliminate r
NumberOfSimulations <- ncol(Simulations)
TargetNumber <- (1-Alpha)*NumberOfSimulations
TargetNumber <- (1 - Alpha) * NumberOfSimulations
KeptSimulations <- Simulations
PresentMinValues <- apply(Simulations, 1, min, na.rm=TRUE)
PresentMaxValues <- apply(Simulations, 1, max, na.rm=TRUE)
PresentMinValues <- apply(Simulations, 1, FUN = min, na.rm = TRUE)
PresentMaxValues <- apply(Simulations, 1, FUN = max, na.rm = TRUE)
# Loop until the target number of simulations is kept
while(ncol(KeptSimulations) > TargetNumber) {
# Remember previous min and max
PreviousMinValues <- PresentMinValues
PreviousMaxValues <- PresentMaxValues
# Select the simulations that gave extreme values
SimulationsToDrop <- c(unlist(apply(KeptSimulations, 1, which.min)), unlist(apply(KeptSimulations, 1, which.max)))
SimulationsToDrop <- c(
unlist(apply(KeptSimulations, MARGIN = 1, FUN = which.min)),
unlist(apply(KeptSimulations, MARGIN = 1, FUN = which.max))
)
# Drop them
KeptSimulations <- KeptSimulations[, -SimulationsToDrop]
# Fails if no simulations are left
if (is.null(dim(KeptSimulations)))
if (is.null(dim(KeptSimulations))) {
stop("Global envelope could not be calculated. More simulations are necessary.")
}
# Calculate min and max
PresentMinValues <- apply(KeptSimulations, 1, min, na.rm=TRUE)
PresentMaxValues <- apply(KeptSimulations, 1, max, na.rm=TRUE)
PresentMinValues <- apply(KeptSimulations, 1, FUN = min, na.rm = TRUE)
PresentMaxValues <- apply(KeptSimulations, 1, FUN = max, na.rm = TRUE)
}
# Interpolate because the kept number of simulations is not always the target
NumberOfKeptSimulations <- ncol(KeptSimulations)
Glo <- PresentMinValues + (PreviousMinValues-PresentMinValues)/NumberOfSimulations*(TargetNumber-NumberOfKeptSimulations)
Ghi <- PresentMaxValues + (PreviousMaxValues-PresentMaxValues)/NumberOfSimulations*(TargetNumber-NumberOfKeptSimulations)
Glo <- PresentMinValues +
(PreviousMinValues - PresentMinValues) / NumberOfSimulations *
(TargetNumber - NumberOfKeptSimulations)
Ghi <- PresentMaxValues +
(PreviousMaxValues - PresentMaxValues) / NumberOfSimulations *
(TargetNumber - NumberOfKeptSimulations)
return(rbind(Glo, Ghi))
}
29 changes: 18 additions & 11 deletions R/GoFtest.R
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
GoFtest <-
function(Envelope) {
GoFtest <- function(Envelope) {

# Verify Envelope
if (!inherits(Envelope, "envelope"))
if (!inherits(Envelope, "envelope")) {
stop("Envelope is not of class envelope")
}
# Verify simulations
if (is.null(attr(Envelope, "simfuns"))) {
stop("Envelope does not contain simulations in its attribute simfuns")
Expand All @@ -14,24 +14,31 @@ function(Envelope) {
}

NumberOfSimulations <- dim(SimulatedValues)[2]
AverageSimulatedValues <- apply(SimulatedValues, 1, sum)/(NumberOfSimulations-1)
rIncrements <- (r-c(0,r)[seq_along(r)])[-1]
AverageSimulatedValues <- apply(SimulatedValues, 1, sum) / (NumberOfSimulations - 1)
rIncrements <- (r - c(0,r)[seq_along(r)])[-1]

# Ui calculate the statistic for a simulation
Ui <- function(SimulationNumber) {
Departure <- (SimulatedValues[, SimulationNumber]-AverageSimulatedValues)[seq_along(r)-1]
WeightedDeparture <- (Departure[!is.nan(Departure)])^2*rIncrements[!is.nan(Departure)]
Departure <- (SimulatedValues[, SimulationNumber] -
AverageSimulatedValues)[seq_along(r) - 1]
WeightedDeparture <- (Departure[!is.nan(Departure)])^2 *
rIncrements[!is.nan(Departure)]
return(sum(WeightedDeparture))
}

# Calculate the Ui statistic for all simulations
SimulatedU <- vapply(seq_len(NumberOfSimulations), Ui, FUN.VALUE=0.0)
SimulatedU <- vapply(
seq_len(NumberOfSimulations),
FUN = Ui,
FUN.VALUE = 0
)

# Calculate the statistic for the actual value
RecenteredValues <- (ActualValues-AverageSimulatedValues)[seq_along(r)-1]
WeightedRecenteredValues <- (RecenteredValues[!is.nan(RecenteredValues)])^2*rIncrements[!is.nan(RecenteredValues)]
RecenteredValues <- (ActualValues - AverageSimulatedValues)[seq_along(r) - 1]
WeightedRecenteredValues <- (RecenteredValues[!is.nan(RecenteredValues)])^2 *
rIncrements[!is.nan(RecenteredValues)]
ActualU <- sum(WeightedRecenteredValues)

# Return the rank
return(mean(ActualU<SimulatedU))
return(mean(ActualU < SimulatedU))
}
70 changes: 49 additions & 21 deletions R/KEnvelope.R
Original file line number Diff line number Diff line change
@@ -1,31 +1,59 @@
KEnvelope <-
function(X, r = NULL, NumberOfSimulations = 100, Alpha = 0.05,
ReferenceType = "", NeighborType = ReferenceType, SimulationType = "RandomPosition",
Precision = 0, Global = FALSE, verbose = interactive()) {
KEnvelope <- function(
X,
r = NULL,
NumberOfSimulations = 100,
Alpha = 0.05,
ReferenceType = "",
NeighborType = ReferenceType,
SimulationType = "RandomPosition",
Precision = 0,
Global = FALSE,
verbose = interactive()) {

CheckdbmssArguments()

# Choose the null hypothesis
SimulatedPP <- switch (SimulationType,
RandomPosition = expression(rRandomPositionK(X, Precision=Precision, CheckArguments = FALSE)),
RandomLabeling = expression(rRandomLabeling(X, CheckArguments = FALSE)),
PopulationIndependence = expression(rPopulationIndependenceK(X, ReferenceType, NeighborType, CheckArguments = FALSE))
)
if (is.null(SimulatedPP))
SimulatedPP <- switch (
SimulationType,
RandomPosition = expression(
rRandomPositionK(X, Precision = Precision, CheckArguments = FALSE)
),
RandomLabeling = expression(
rRandomLabeling(X, CheckArguments = FALSE)
),
PopulationIndependence = expression(
rPopulationIndependenceK(
X,
ReferenceType,
NeighborType,
CheckArguments = FALSE)
)
)
if (is.null(SimulatedPP)) {
stop(paste("The null hypothesis", sQuote(SimulationType), "has not been recognized."))
}
# local envelope, keep extreme values for lo and hi (nrank=1)
Envelope <- envelope(X, fun=Khat, nsim=NumberOfSimulations, nrank=1,
r=r, ReferenceType=ReferenceType, NeighborType=NeighborType,
CheckArguments = FALSE,
simulate=SimulatedPP, verbose=verbose, savefuns=TRUE
)
attr(Envelope, "einfo")$H0 <- switch (SimulationType,
RandomPosition = "Random Position",
RandomLabeling = "Random Labeling",
PopulationIndependence = "Population Independence"
)
Envelope <- envelope(
X,
fun = Khat,
nsim = NumberOfSimulations,
nrank = 1,
r = r,
ReferenceType = ReferenceType,
NeighborType = NeighborType,
CheckArguments = FALSE,
simulate = SimulatedPP,
verbose = verbose,
savefuns = TRUE
)
attr(Envelope, "einfo")$H0 <- switch (
SimulationType,
RandomPosition = "Random Position",
RandomLabeling = "Random Labeling",
PopulationIndependence = "Population Independence"
)
# Calculate confidence intervals
Envelope <- FillEnvelope(Envelope, Alpha, Global)
Envelope <- FillEnvelope(Envelope, Alpha = Alpha, Global = Global)
# Return the envelope
return (Envelope)
}
77 changes: 53 additions & 24 deletions R/KdEnvelope.R
Original file line number Diff line number Diff line change
@@ -1,34 +1,63 @@
KdEnvelope <-
function(X, r = NULL, NumberOfSimulations = 100, Alpha = 0.05,
ReferenceType, NeighborType = ReferenceType, Weighted = FALSE, Original = TRUE,
Approximate = ifelse(X$n < 10000, 0, 1), Adjust = 1, MaxRange = "ThirdW",
StartFromMinR = FALSE,
SimulationType = "RandomLocation", Global = FALSE, verbose = interactive()) {
KdEnvelope <- function(
X,
r = NULL,
NumberOfSimulations = 100,
Alpha = 0.05,
ReferenceType,
NeighborType = ReferenceType,
Weighted = FALSE,
Original = TRUE,
Approximate = ifelse(X$n < 10000, 0, 1),
Adjust = 1,
MaxRange = "ThirdW",
StartFromMinR = FALSE,
SimulationType = "RandomLocation",
Global = FALSE,
verbose = interactive()) {

CheckdbmssArguments()

# Choose the null hypothesis
SimulatedPP <- switch (SimulationType,
RandomLocation = expression(rRandomLocation(X, CheckArguments = FALSE)),
RandomLabeling = expression(rRandomLabelingM(X, CheckArguments = FALSE)),
PopulationIndependence = expression(rPopulationIndependenceM(X, ReferenceType, CheckArguments = FALSE))
)
if (is.null(SimulatedPP))
SimulatedPP <- switch (
SimulationType,
RandomLocation = expression(rRandomLocation(X, CheckArguments = FALSE)),
RandomLabeling = expression(rRandomLabelingM(X, CheckArguments = FALSE)),
PopulationIndependence = expression(
rPopulationIndependenceM(X, ReferenceType, CheckArguments = FALSE)
)
)
if (is.null(SimulatedPP)) {
stop(paste("The null hypothesis", sQuote(SimulationType), "has not been recognized."))
}

# local envelope, keep extreme values for lo and hi (nrank=1)
Envelope <- envelope(X, fun=Kdhat, nsim=NumberOfSimulations, nrank=1,
r=r, ReferenceType=ReferenceType, NeighborType=NeighborType, Weighted=Weighted, Original=Original,
Approximate=Approximate, Adjust=Adjust, MaxRange=MaxRange, StartFromMinR=StartFromMinR,
CheckArguments = FALSE,
simulate=SimulatedPP, verbose=verbose, savefuns=TRUE
)
attr(Envelope, "einfo")$H0 <- switch (SimulationType,
RandomLocation = "Random Location",
RandomLabeling = "Random Labeling",
PopulationIndependence = "Population Independence"
)
Envelope <- envelope(
X,
fun = Kdhat,
nsim = NumberOfSimulations,
nrank = 1,
r = r,
ReferenceType = ReferenceType,
NeighborType = NeighborType,
Weighted = Weighted,
Original = Original,
Approximate = Approximate,
Adjust = Adjust,
MaxRange = MaxRange,
StartFromMinR = StartFromMinR,
CheckArguments = FALSE,
simulate = SimulatedPP,
verbose = verbose,
savefuns = TRUE
)
attr(Envelope, "einfo")$H0 <- switch (
SimulationType,
RandomLocation = "Random Location",
RandomLabeling = "Random Labeling",
PopulationIndependence = "Population Independence"
)
# Calculate confidence intervals
Envelope <- FillEnvelope(Envelope, Alpha, Global)
Envelope <- FillEnvelope(Envelope, Alpha = Alpha, Global = Global)
# No edge effect correction
attr(Envelope, "einfo")$valname <- NULL
# Return the envelope
Expand Down
24 changes: 18 additions & 6 deletions R/Khat.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
Khat <-
function(X, r = NULL, ReferenceType = "", NeighborType = ReferenceType, CheckArguments = TRUE) {
Khat <- function(
X,
r = NULL,
ReferenceType = "",
NeighborType = ReferenceType,
CheckArguments = TRUE) {

if (CheckArguments) {
CheckdbmssArguments()
Expand All @@ -11,17 +15,25 @@ function(X, r = NULL, ReferenceType = "", NeighborType = ReferenceType, CheckArg

# K intra
if (ReferenceType == "" & NeighborType == "") {
return (Kest(X, r=r, correction="best"))
return (Kest(X, r = r, correction = "best"))
}
# K intra for a single point type
if (ReferenceType == NeighborType) {
X.reduced <- X[marks(X)$PointType == ReferenceType]
return (Kest(X.reduced, r=r, correction="best"))
return (Kest(X.reduced, r = r, correction = "best"))
}
# K inter calls Kcross. The marks must contain the type, with no weight.
if (ReferenceType != NeighborType) {
X.cross <- X
spatstat.geom::marks(X.cross) <- marks(X)$PointType
return (Kcross(X.cross, i=ReferenceType , j=NeighborType, r=r, correction="best"))
marks(X.cross) <- marks(X)$PointType
return(
Kcross(
X.cross,
i = ReferenceType ,
j = NeighborType,
r = r,
correction = "best"
)
)
}
}
Loading

0 comments on commit d4f3504

Please sign in to comment.