Closest-pair problem: Difference between revisions

Using dist function for brute force, but divide and conquer (as per pseudocode) for speed:

<lang 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)

}

xp <- x[order(x[,"x"]),]

yp <- x[order(x[,"y"]),]

.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,]

xm <- xp[floor(n/2),"x"]

yl <- yp[which(yp[,"x"] <= xm),]

yr <- yp[which(yp[,"x"] > xm),]

cpl <- .cpdandc.rec(xl,yl)

cpr <- .cpdandc.rec(xr,yr)

if (cpl$d<cpr$d) cp <- cpl else cp <- cpr

dmin <- cp$d

ys <- yp[which(abs(xm - yp[,"x"]) < dmin),]

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"] < dmin)

{

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

}

}

}

cp

}

}

.cpdandc.rec(xp,yp)

}

# 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))

</lang>

Sample output

<pre>

How many points?

1: 500

Read 1 item

Closest pairs divide and conquer:

$p1

x y

1.68807938 0.05876328

$p2

x y

1.68904694 0.05878173

$d

[1] 0.0009677302

That took 0.43 seconds.

Closest pairs brute force:

$p1

x y

1.68807938 0.05876328

$p2

x y

1.68904694 0.05878173

$d

[1] 0.0009677302

That took 6.38 seconds.

</pre>