
N <- 10000
D <- 1
coef.1 <- matrix(NA,N,2)
coef.2 <- matrix(NA,N,2)
path.1 <- matrix(NA,N,2)
path.2 <- matrix(NA,N,2)
path.1[1,] <- c(40,40)
path.2[1,] <- c(60,60)
d.start <- sqrt(sum((path.1[1,]-path.2[1,])^2))
ch <- "."
set.seed(1)
system.time({
for (i in 2:N){
if (i%%50==0) cat(ch)
path.1[i,] <- path.1[i-1,] + sample(-D:D,2)
path.2[i,] <- path.2[i-1,] + sample(-D:D,2)
coef.1[i,] <- get.line(path.1[(i-1):i,])
coef.2[i,] <- get.line(path.2[(i-1):i,])
r.1 <- sqrt(max(rowSums((path.1[1:i,]-path.1[1,])^2)))
r.2 <- sqrt(max(rowSums((path.2[1:i,]-path.2[1,])^2)))
if (r.1+r.2 < d.start) next # paths do not overlap
ch <- "1"
d.1 <- sqrt(min(rowSums((path.2[1:i,]-path.1[1,])^2)))
d.2 <- sqrt(min(rowSums((path.1[1:i,]-path.2[1,])^2)))
if (d.1>r.1 & d.2>r.2) next
ch <- "2"
cross <- sapply(2:i,
function(k){seg.intersect(path.2[(k-1):k,],path.1[(i-1):i,],k)})
if (any(cross)) break
cross <- sapply(2:i,
function(k){seg.intersect(path.1[(k-1):k,],path.2[(i-1):i,],k)})
if (any(cross)) break
}
})
# 11111111112222222222222222222222
# user system elapsed
# 1016.82 0.13 1024.18
print(paste("End at Step: ",i))
# [1] "End at Step: 1624"
plot(0:100,0:100, type="n", xlab="X", ylab="Y")
points(path.1[1,1],path.1[1,2], pch=16, col="red", cex=1.5)
points(path.2[1,1],path.2[1,2], pch=16, col="green", cex=1.5)
lines(path.1[1:i,])
lines(path.2[1:i,],col="red")
As @CarlWitthoft points out, at each step you must check all the previous line segments for crossings. This creates a serious problem, because at each new step, i
, there are 2*(i-1)
tests for crossings. So, if you get to a crossing at step k
, there will have been 2*k*(k+1)
tests. If k ~O(10000)
, then there can be potentially 100MM tests.
To make this more efficient, we store not only the two new points at each step, but also the slope and intercept of the newly created line segments. This avoids recalculating slope and intercept for all prior line segments at each step. In addition, we calculate the path radius, r, for each path at each step. This is the distance between the starting point and the point on the path farthest from the starting point. If the distance between the path starting point and the closest point on the other path is larger than the path radius, there can be no crossings and we can skip the individual segment comparisons for this step.
Your problem is interesting for other reasons. The normal way to test for crossings is to determine if the intersection between the two lines is on either segment. This is cumbersome but straightforward. However there are many special cases: Are the lines parallel? If so, are they coincident? If so do the segments overlap? What about vertical lines (slope=Inf). Because you set the increment to a random integer on [-1,1], all of these possibilities are quite likely to happen eventually in a path with 10000 steps. So the function seg.intersect(...)
above has to account for all these possibilities. You would think there is a function in R that does that, but I couldn't find one, so here is a (messy) version:
get.line <- function(l) { # returns slope and intercept
if (diff(l)[1]==0) return(c(Inf,NA))
m <- diff(l)[2]/diff(l)[1]
b <- l[1,2]-m*l[1,1]
return(c(m,b))
}
is.between <- function(x,vec) { # test if x is between values in vec
return(x>=range(vec)[1] & x<=range(vec)[2])
}
special.cases = function(l1,l2, coeff) {
# points coincide: not a line segment!
if (rowSums(diff(l1)^2)==0 | rowSums(diff(l2)^2)==0) return(c(NA,FALSE))
# both lines vertical
if (is.infinite(coeff[1,1]) & is.infinite(coeff[2,1])) {
if (l1[1,1]!=l2[1,1]) return(c(NA,FALSE))
t1 <- is.between(l1[1,2],l2[,2]) | is.between(l1[2,2],l2[,2])
t2 <- is.between(l2[1,2],l1[,2]) | is.between(l2[2,2],l1[,2])
return(c(NA,t1|t2))
}
# only l1 is vertical
if (is.infinite(coeff[1,1]) & is.finite(coeff[2,1])) {
x <- l1[1,1]
y <- c(x,1) %*% coeff[2,]
return(c(x,y))
}
# only l2 is vertical
if (is.finite(coeff[1,1]) & is.infinite(coeff[2,1])) {
x <- l2[1,1]
y <- c(x,1) %*% coeff[1,]
return(c(x,y))
}
# parallel, non-coincident lines
if (diff(coeff[,1])==0 & diff(coeff[,2])!=0) return(c(NA,FALSE))
# parallel, coincident lines
if (diff(coeff[,1])==0 & diff(coeff[,2])==0) {
x <- l1[1,1]
y <- l1[1,2]
return(c(x,y))
}
# base case: finite slopes, not parallel
x <- -diff(coeff[,2])/diff(coeff[,1])
y <- c(x,1) %*% coeff[1,]
return(c(x,y))
}
seg.intersect <- function(l1,l2,i){
pts <- list(l1,l2)
coeff <- rbind(coef.1[i,],coef.2[i,])
z <- special.cases(l1,l2, coeff)
if (is.na(z[1])) return (z[2])
# print(coeff)
# print(z)
found <- do.call("&",
lapply(pts,function(x){is.between(z[1],x[,1]) & is.between(z[2],x[,2])}))
return(found)
}