88 lines
1.7 KiB
R
88 lines
1.7 KiB
R
closest.pairs.bruteforce <- function(x, y=NULL)
|
|
{
|
|
if (!is.null(y))
|
|
{
|
|
x <- cbind(x,y)
|
|
}
|
|
d <- dist(x)
|
|
cp <- x[combn(1:nrow(x), 2)[, which.min(d)],]
|
|
list(p1=cp[1,], p2=cp[2,], d=min(d))
|
|
}
|
|
|
|
closest.pairs.dandc <- function(x, y=NULL)
|
|
{
|
|
if (!is.null(y))
|
|
{
|
|
x <- cbind(x,y)
|
|
}
|
|
if (sd(x[,"x"]) < sd(x[,"y"]))
|
|
{
|
|
x <- cbind(x=x[,"y"],y=x[,"x"])
|
|
swap <- TRUE
|
|
}
|
|
else
|
|
{
|
|
swap <- FALSE
|
|
}
|
|
xp <- x[order(x[,"x"]),]
|
|
.cpdandc.rec <- function(xp,yp)
|
|
{
|
|
n <- dim(xp)[1]
|
|
if (n <= 4)
|
|
{
|
|
closest.pairs.bruteforce(xp)
|
|
}
|
|
else
|
|
{
|
|
xl <- xp[1:floor(n/2),]
|
|
xr <- xp[(floor(n/2)+1):n,]
|
|
cpl <- .cpdandc.rec(xl)
|
|
cpr <- .cpdandc.rec(xr)
|
|
if (cpl$d<cpr$d) cp <- cpl else cp <- cpr
|
|
cp
|
|
}
|
|
}
|
|
cp <- .cpdandc.rec(xp)
|
|
|
|
yp <- x[order(x[,"y"]),]
|
|
xm <- xp[floor(dim(xp)[1]/2),"x"]
|
|
ys <- yp[which(abs(xm - yp[,"x"]) <= cp$d),]
|
|
nys <- dim(ys)[1]
|
|
if (!is.null(nys) && nys > 1)
|
|
{
|
|
for (i in 1:(nys-1))
|
|
{
|
|
k <- i + 1
|
|
while (k <= nys && ys[i,"y"] - ys[k,"y"] < cp$d)
|
|
{
|
|
d <- sqrt((ys[k,"x"]-ys[i,"x"])^2 + (ys[k,"y"]-ys[i,"y"])^2)
|
|
if (d < cp$d) cp <- list(p1=ys[i,],p2=ys[k,],d=d)
|
|
k <- k + 1
|
|
}
|
|
}
|
|
}
|
|
if (swap)
|
|
{
|
|
list(p1=cbind(x=cp$p1["y"],y=cp$p1["x"]),p2=cbind(x=cp$p2["y"],y=cp$p2["x"]),d=cp$d)
|
|
}
|
|
else
|
|
{
|
|
cp
|
|
}
|
|
}
|
|
|
|
# Test functions
|
|
cat("How many points?\n")
|
|
n <- scan(what=integer(),n=1)
|
|
x <- rnorm(n)
|
|
y <- rnorm(n)
|
|
tstart <- proc.time()[3]
|
|
cat("Closest pairs divide and conquer:\n")
|
|
print(cp <- closest.pairs.dandc(x,y))
|
|
cat(sprintf("That took %.2f seconds.\n",proc.time()[3] - tstart))
|
|
plot(x,y)
|
|
points(c(cp$p1["x"],cp$p2["x"]),c(cp$p1["y"],cp$p2["y"]),col="red")
|
|
tstart <- proc.time()[3]
|
|
cat("\nClosest pairs brute force:\n")
|
|
print(closest.pairs.bruteforce(x,y))
|
|
cat(sprintf("That took %.2f seconds.\n",proc.time()[3] - tstart))
|