2008-05-18

線分データの数え上げ(続きの続きの続き)

計算幾何学探索木 R

前の記事はこちら
今、N個のデータが得られるとする。そのデータは、１次元尺度上の線分に相当するデータであるものとする。簡単に言えば、aiからbiまでというようなデータが、N個得られるとする。このようなデータから、ある値xについて、[tex:ai<=c

test<-function(file){
d<-read.table(file=file)
start<-c(d$V1)
end<-c(d$V2)
startend<-c(start,end)
value<-c(rep(1,length(start)),rep(-1,length(end)))
mat<-matrix(c(startend,value),ncol=2)
matb<-matrix(rep(1,length(startend)^2),ncol=length(startend))
matb[upper.tri(matb)]<-0
order<-order(mat[,1])
accum<-matb %*% mat[order,2]
retmat<-matrix(c(mat[order,1],accum),ncol=2)
return(retmat)
}

行列を使わずにループを使うと

test2<-function(file){
d<-read.table(file=file)
start<-c(d$V1)
end<-c(d$V2)
startend<-c(start,end)
value<-c(rep(1,length(start)),rep(-1,length(end)))
mat<-matrix(c(startend,value),ncol=2)
order<-order(mat[,1])
accum<-rep(1,length(startend))
accum[1]=mat[order,2][1]
for(i in 2:length(startend)){
accum[i]<-accum[i-1]+mat[order,2][i]
}

retmat<-matrix(c(mat[order,1],accum),ncol=2)
return(retmat)
}

線分の始点・終点の向きを考慮しない場合と、考慮してある向きだけについて加算する、などとするための条件を入れると

test3signPN<-function(file,pn){
d<-read.table(file=file)
start<-c(d$V1)
end<-c(d$V2)
if(pn<0){
sign<-sign(start-end)
}else if(pn>0){
sign<-sign(end-start)
}else{
sign<-rep(1,length(start))
}
sign2<-c(sign,sign)
sign2<-(sign2+1)/2
startend<-c(start,end)
value<-c(rep(1,length(start)),rep(-1,length(end)))
value<-value*sign2
mat<-matrix(c(startend,value),ncol=2)
order<-order(mat[,1])
accum<-rep(1,length(startend))
accum[1]=mat[order,2][1]
for(i in 2:length(startend)){
accum[i]<-accum[i-1]+mat[order,2][i]
}

retmat<-matrix(c(mat[order,1],accum),ncol=2)
return(retmat)
}

データがテキストファイル "test.txt"に、１データ１行のタブ区切りで次のように入っているとする

Rでこれを読んで処理するとすると・・・

#ファイルを読む
> d<-read.table(file="test.txt")
#読んだファイルを見てみる
> d
  V1 V2
1  4  6
2  2  5
3  7  9
4  1  8
#始点はV1のカラムなので、それを取り出す
> start<-c(d$V1)
> start
[1] 4 2 7 1
#終点はV2のカラムなので、それを取り出す
> end<-c(d$V2)
#始点と終点を同じベクトルに入れる
> startend<-c(start,end)
> startend
[1] 4 2 7 1 6 5 9 8
#始点と終点を併せた要素数と同じ長さで、始点には1が終点には-1が対応するような、ベクトルを作る
> value<-c(rep(1,length(start)),rep(-1,length(end)))
> value
[1]  1  1  1  1 -1 -1 -1 -1
#始点終点の位置を第一列に持ち、それに与えた1,-1の値を第２列に持つような行列を作る
> mat<-matrix(c(startend,value),ncol=2)
> mat
     [,1] [,2]
[1,]    4    1
[2,]    2    1
[3,]    7    1
[4,]    1    1
[5,]    6   -1
[6,]    5   -1
[7,]    9   -1
[8,]    8   -1
#始点終点の位置でソートする
> mat[order(mat[,1]),]
     [,1] [,2]
[1,]    1    1
[2,]    2    1
[3,]    4    1
[4,]    5   -1
[5,]    6   -1
[6,]    7    1
[7,]    8   -1
[8,]    9   -1
#始点終点の位置でソートしたときの第２列ベクトルを確認する
> mat[order(mat[,1]),2]
[1]  1  1  1 -1 -1  1 -1 -1
#インクリメント・デクリメント作業のために、三角行列を作る
> matb<-matrix(rep(1,length(startend)^2),ncol=length(startend))
> matb
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    1    1    1    1    1    1    1
[2,]    1    1    1    1    1    1    1    1
[3,]    1    1    1    1    1    1    1    1
[4,]    1    1    1    1    1    1    1    1
[5,]    1    1    1    1    1    1    1    1
[6,]    1    1    1    1    1    1    1    1
[7,]    1    1    1    1    1    1    1    1
[8,]    1    1    1    1    1    1    1    1
> matb[upper.tri(matb)]<-0
> matb
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    0    0    0    0    0
[2,]    1    1    0    0    0    0    0    0
[3,]    1    1    1    0    0    0    0    0
[4,]    1    1    1    1    0    0    0    0
[5,]    1    1    1    1    1    0    0    0
[6,]    1    1    1    1    1    1    0    0
[7,]    1    1    1    1    1    1    1    0
[8,]    1    1    1    1    1    1    1    1
#インクリメント・デクリメントの結果を行列の積として計算する
> matb %*% mat[order(mat[,1]),2]
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    2
[5,]    1
[6,]    2
[7,]    1
[8,]    0