十字路口的两步功能 - intersection of two step functions

- 此内容更新于:2016-02-24
主题:

我试图确定点(x,y)两个函数相交的地方。之间的步插值函数的点集。一个函数是弱增加()和其他弱减少()。我R中的编码,但一般的算法也是好的。如果有帮助,这是确定市场供需平衡集点。两个向量的长度是不同的,他们的x和y不会是相同的。一些示例数据:在这个例子中,十字路口在(x=2.7275363,y=2.7275363),来自v2和从v1。谢谢

原文:

I trying to determine the point (x,y) where two functions intersect. The functions are the step interpolation between sets of points. One function is weakly increasing (v1) and the other weakly decreasing (v2). I'm coding in R, but a general algorithm is also ok.

If it helps, this is to determine market equilibrium with sets of supply and demand points.

The length of the two vectors is different and their x's and y's will not be the same.

Some example data:

set.seed(4)

v1 = data.frame( y = cumsum( runif(10) ) ,
                 x = cumsum( runif(10) ) )
v2 = data.frame( y = 5-cumsum( runif(8) )  ,
                 x = cumsum( runif(8) ) )

plot(y=0,x=0,type="n",xlim=c(0,5),ylim=c(0,5),xlab="x",ylab="y")

lines( y=v1$y , x=v1$x , type="S" , col="blue" )
lines( y=v1$y , x=v1$x , type="p" , col="blue" )

lines( y=v2$y , x=v2$x , type="s" , col="red" )
lines( y=v2$y , x=v2$x , type="p" , col="red" )

In this example, the intersection is at (x=2.7275363 , y=2.510405), where the x is from v2 and y is from v1.

Thanks

网友:其实确切的交点似乎在(x=2.727536,y=2.727536)。你可以确认你希望的答案是什么?

(原文:Actually the exact point of intersection seem to to be at (x=2.727536 , y=2.239863). You can confirm what you expect the answer to be?)

楼主:我离开我的电脑在工作所以不能证实,但它可能是我犯了一个错误的y号码。所以,把它作为。我的歉意。

(原文:I've left my computer at work so can't confirm, but it is likely I made a mistake with the y number. So, take it as being y=2.239863. My apologies.)

解决方案:
你画线不同在每种情况下步:v1首先你改变垂直,然后水平(和),而对于v2你反向顺序(跨)。假设这是正确的,那么你的交点后将立即或在v1下点沿轴是v1y坐标较低。我们可以发现通过:现在有两个不同的案例——十字路口可以水平或垂直臂下面这一点从v1。它将立即前如果事先v2高于v1后我们发现,否则它将在水平的手臂。这是显而易见的如果你画出来,我会附上图片如果你没有看到这一点。令人担忧的是,这并不同意你(x=2.7275363,y=2.510405)的答案,但当我的阴谋,我出现在十字路口。所以:我没有明白你想要什么;你错误,或者有一个不同的方案对水平和垂直组件的顺序。上面的代码应该适应不同的方案。
原文:

You're drawing your step lines differently in each case: v1 you change the vertical first, and then the horizontal (up and across), whereas for v2 you reverse the order (across then down). Assuming this is correct, then your intersection point will be at or immediately after a point in v1 where the next point along the axis is a v1 with a lower y coordinate. We can find that by doing:

v1$v <- 1
v2$v <- 2
v3 <- rbind(v1,v2)
v3 <- v3[order(v3$x),]
v3$diff <- c( diff(v3$y),0)
ind <- which(v3$diff < 0 & v3$v ==1)[1]

There are now two distinct cases - the intersection could be on the horizontal or vertical arm following this point from v1. It will be the former if the immediately preceeding v2 is higher than the v1 after our found one; otherwise it will be in the horizontal arm. This is clear if you draw it out - I'll try and attach an image if you don't see this.

previousV2 <- tail(which(v3$v[1:ind]==2),1)
nextV1 <- which(v3$v[-(1:ind)]==1)[1] + ind
if (v3$y[previousV2] > v3$y[nextV1]) {
  x <- v3$x[ind+1]
  y <- v3$y[nextV1]
} else {
  x <- v3$x[ind]
  y <- v3$y[previousV2]
}

Worryingly, this doesn't agree with your (x=2.7275363 , y=2.510405) answer, but when I plot it, mine appears on the intersection. So either: I haven't understood what you want; you've miscalculated; or there's a different scheme regarding the order of horizontal and vertical components. The above code should be adaptable to different schemes.

网友:尝试你的方法和运行。我添加到最终结果。似乎事情可能会变得很棘手。必须与R决定如何画线。如果关键是找到的交集R画的线,这是一回事,但是也许这不是答案OP之后。希望他可以澄清。

(原文:Try running your method with set.seed(5) and set.seed(7). I added abline(v=x, lty=2); abline(h=y, lty=2) to the end to see the results. It seems like things can get tricky. Must have to do with how R decides to draw the lines. If the point is to find the intersection of the lines that R draws, that's one thing; but maybe that's not really the answer the OP is after. Hopefully he can clarify.)

网友:嗯,好的。应该有一些逻辑抓的这种行为。我认为我们需要一个,看到印第安纳州发生的第一次——是一个漫长的一周,和不能快速看看这个条件需要它自己的逻辑所示我的第二块的代码,但是希望OP,它应该提供一个依据

(原文:Hmmm, good catch. There should be some logic that catches this behaviour. I think we need an ind2 <- which(v3$diff > 0 & v3$v==2) and see which of the ind's occurs first- been a long week, and can't quickly see if this conditional needs it's own logic shown in my second chunk of code or not, but hopefully it should provide a basis for the OP)

网友:关于R如何绘制线条,这取决于参数的情况下是否水平或垂直第一——我认为这是一个深思熟虑的区别这个问题,但也许不是。

(原文:Regarding how R draws the lines, it depends on the case of the type="s/S" parameter whether it's horizontal or vertical first - I assumed this was a deliberate distinction in the question, but maybe not.)

网友:我错过了“s/s”的区别可能是故意的。眼睛看起来很简单,但是有点欺骗者的代码。但是你比我有提供一个更干净的开始。

(原文:I missed the "s/S" distinction so likely that was intentional. It seems so simple by eye, but is a bit tricker by code. But you did provide a cleaner start than I had.)

楼主:是的,s/s的区别是故意的。

(原文:Yes, the s/S distinction was deliberate.)

解决方案:
我似乎有东西,但它是一个更复杂的比我期望。首先,让我定义一个helper函数这只是帮助检查是否两人之间的一个数字。现在我将这两个数据。框架集的连续点对,然后我检查每一个可能的组合段重叠在正确的方式。(这很重要,这是“S”的)。现在包含匹配点集,我只需要提取正确的坐标取决于哪种类型的十字路口。也许有一个更好的办法,但至少你有另一个方法似乎比较答案与工作。
原文:

I seem to have something that works but it's a lot more complicated than i was expecting.

First, let me define a helper function

between <- function(x, a, b) {
    if(missing(b)) {
        if(length(a)==2) {
          a<-t(a)
       }
    } else {
        a <- unname(cbind(a,b))
    }
    a<-t(apply(a,1,sort))
    a[,1] <= x & x <= a[,2]
}

this just helps to check if a number is between two others. Now I will embed the two data.frames to make sets of consecutive point pairs, then i check each possible combination for segments that overlap in just the right way. (It's important that v1 here is the "S" and v2 is the s.)

sa<-embed(as.matrix(v1[,c("x","y")]),2)
sz<-embed(as.matrix(v2[,c("x","y")]),2)
xx<-outer(1:nrow(sa), 1:nrow(sz), function(a,z)
    (between(sa[a,2], sz[z,c(2,4)]) & between(sz[z,1], sa[a,c(1,3)])) *1 
     + (between(sz[z,4], sa[a,c(2,4)]) & between(sa[a,3], sz[z,c(1,3)]))*2
)

Now xx contains the matching set of points, I just need to extract the correct coordinates depending on which type of intersection occurred.

i <- which(xx!=0, arr.ind=T)
int.pt <- if(nrow(i)>0 && ncol(i)==2) {
    if(xx[i]==1) {
       c(sz[i[2],1], sa[i[1],2])
    } else if (xx[i]==2) {
       c(sa[i[1],3], sz[i[2],4])
    }
} else {
   c(NA,NA)
}
#optionally plot intersection
#if (all(!is.na(int.pt))) {
#    points(int.pt[1],int.pt[2], pch=20, col="black")
#    abline(v=int.pt[1], h=int.pt[2], lty=2)
#}

Perhaps there is a better way, but at least you have another method that seems to work to compare answers with.

sample intersection detection

楼主:谢谢!这种方法适用。我会观看它,如果我能加速我的解决方案。

(原文:Thanks! This approach works well. I'll play around with it and post my solution if I can speed it up.)

解决方案:
我有另一个思考问题。一个关键问题是,我需要找到交集在一个优化程序,所以它必须快。所以,我想出了以下(包括在未来的其他人有同样的问题)。这是一个修改Bentley-Ottmann算法。
原文:

I had another think about the problem. A key issue is that I need to find the intersection within an optimisation routine, so it has to be fast. So, I came up with the following (included here in case others have to same problem in the future). It is a modified Bentley-Ottmann algorithm.

# create some data
supply = data.frame( p =                    cumsum( runif(1000) ) ,
                     q =                    cumsum( runif(1000) ) )
demand = data.frame( p = tail(supply,1)$p - cumsum( runif(1000) )  ,
                     q =                    cumsum( runif(1000) ) )

# create tables that identify coordinates of horizontal and vertical lines 
demand.h = cbind( p       = head(demand,-1)$p , 
                  q.lower = head(demand,-1)$q , 
                  q.upper = tail(demand,-1)$q )

supply.v = cbind( q       = head(supply,-1)$q ,
                  p.lower = head(supply,-1)$p ,
                  p.upper = tail(supply,-1)$p )

demand.v = cbind( q       = tail(demand,-1)$q ,
                  p.lower = tail(demand,-1)$p ,
                  p.upper = head(demand,-1)$p )

supply.h = cbind( p       = tail(supply,-1)$p ,
                  q.lower = head(supply,-1)$q ,
                  q.upper = tail(supply,-1)$q )

# define a function
find.intersection = function( f.A , f.B ){
  f.result = any( f.B[,2]<=f.A[1]  & f.B[,3]>=f.A[1]  & 
                  f.A[2] <=f.B[,1] & f.A[3] >=f.B[,1] )
  return( f.result )
}

# find the intersection
intersection.h = c( demand.h[ apply( demand.h , 
                                     MARGIN=1 , 
                                     FUN=find.intersection , 
                                     supply.v ) , 1 ] ,
                    supply.v[ apply( supply.v , 
                                     MARGIN=1 , 
                                     FUN=find.intersection , 
                                     demand.h ) , 1 ] )

intersection.v = c( supply.h[ apply( supply.h , 
                                     MARGIN=1 , 
                                     FUN=find.intersection , 
                                     demand.v ) , 1 ] ,
                    demand.v[ apply( demand.v , 
                                     MARGIN=1 , 
                                     FUN=find.intersection , 
                                     supply.h ) , 1 ] )

intersection = c( intersection.h , intersection.v )

# (optional) if you want to print the graph and intersection
plot(y=0,x=0,type="n",
     xlim=c(intersection[2]-1,intersection[2]+1),
     ylim=c(intersection[1]-1,intersection[1]+1),
     xlab="q",ylab="p")

lines( y=supply$p , x=supply$q , type="S" , col="black" )
lines( y=supply$p , x=supply$q , type="p" , col="black" )

lines( y=demand$p , x=demand$q , type="s" , col="black" )
lines( y=demand$p , x=demand$q , type="p" , col="black" )

points(intersection[2],intersection[1], pch=20, col="red")
abline( v=intersection[2], h=intersection[1], lty=2 , col="red")