疎行列処理はアドレス重複が苦手

• 昨日の記事では、曲面の変形・共形変換・スピン変換をRでなぞってみているが、その中で、疎行列を扱っている
• 三角化メッシュの辺の両端頂点IDを番地とする行列を疎行列で作り、その要素を辺が含まれる三角形の数だけ操作をしたいのだが
• 何回の操作になるのかが不定なので、少し厄介
• それをループしないで行うべく、rle()関数などを使ってみる
• まず、疎行列・疎ベクトル
• 空の疎ベクトルを作ってみる。第１要素に0を入れて、長さ１０の疎ベクトルを作る。要素を何も入れないと作れないみたいなので

library(Matrix)
v1 <- sparseVector(c(0),i=c(1),length=10)
v1

> v1
sparse vector (nnz/length = 0/10) of class "dsparseVector"
 [1] . . . . . . . . . .

• 疎行列にしてみる

m1 <- Matrix(v1,2,5)
m1

> m1
2 x 5 sparse Matrix of class "dgTMatrix"
              
[1,] . . . . .
[2,] . . . . .

• 番地制御さえできれば、うまくいきそう

v2 <- sparseVector(c(1,2,3),i=c(3,2,5),length=10)
v2

> v2
sparse vector (nnz/length = 3/10) of class "dsparseVector"
 [1] . 2 1 . 3 . . . . .

• 番地の指定に重複はあってはいけない

v3 <- sparseVector(c(1,2,3),i=c(1,3,1),length=10)

• 以下のようなエラーが出る。このエラーは「昇順」になっていない、からではなくて、重複番地があるから。

> v3 <- sparseVector(c(1,2,3),i=c(1,3,1),length=10)
 以下にエラー validObject(.Object) : 
  invalid class “dsparseVector” object: 'i' must be sorted strictly increasingly

• さて、番地ごとに何回かずつ、アクセスして、加算したいとしよう
  • banchiにvalueを加えたい。以下の例だと、第１番地に+1,+4+8、第２番地に+2+5+10とか

value <- 1:10
banchi <- c(1,2,3,1,2,4,50,1,3,2)

• rle()関数
  • rle()関数は、連続値がいくつあるかを教えてくれる関数

v1 <- c(1,1,1,2,2,3,4,4,4)
rle(v1)
v1.sh <- sample(v1)
v1.sh
rle(v1.sh)

• 以下に見られるように、「ソートすれば、『出現回数』」がわかり、そうでなければ、「続き回数」がわかる

> v1 <- c(1,1,1,2,2,3,4,4,4)
> rle(v1)
Run Length Encoding
  lengths: int [1:4] 3 2 1 3
  values : num [1:4] 1 2 3 4
> v1.sh <- sample(v1)
> v1.sh
[1] 4 1 3 4 2 2 1 1 4
> rle(v1.sh)
Run Length Encoding
  lengths: int [1:7] 1 1 1 1 2 2 1
  values : num [1:7] 4 1 3 4 2 1 4

• 今、やりたいのは、複数回出る疎ベクトル番地には１回のアクセスをすることなので、番地をソートしてからやるのが大事
• 番地のソート順を出して、番地と値とをその順序にする
• 番地をソートしてrle()をかければ、それぞれの番地の出現回数が出る
• 番地数ｘ最大出現回数の行列に値を納めてやって、行ごとに加算すれば、番地ごとの加算数が合算できるから、そのような行列を作りたい
• あくまでも「疎ベクトル」「疎行列」でやりたいので、rle()の出力から、「行列」をベクトル扱いした場合の「番地」を作ってやることにする

my.vector.access <- function(v,a,func=sum,zero=0){
	if(is.vector(a)){
		v <- matrix(v,ncol=1)
	}
	ord <- order(a)
	rle.out <- rle(a[ord])
	num.row <- length(rle.out[[1]])
	num.col <- max(rle.out[[1]])
	tmp1 <- rep(1:num.row,rle.out[[1]])
	tmp2 <- c()
	for(i in 1:num.row){
		tmp2 <- c(tmp2,1:rle.out[[1]][i])
	}
	addr <- tmp1 + num.row*(tmp2-1)
	ret.v <- matrix(0,num.row,ncol(v))
	for(i in 1:ncol(v)){
		if(zero==0){
			tmp.V <- sparseVector(v[ord,i],i=addr,length=num.row*num.col)
			M <- Matrix(tmp.V,num.row,num.col)
		}else{
			M <- matrix(zero,num.row,num.col)
			M[addr] <- v[ord,i]

		}
		ret.v[,i] <- apply(M,1,func)

	}
	return(list(A = rle.out[[2]],V = ret.v))
}

v <- 1:10
a <- banchi <- c(1,2,3,1,2,4,50,1,3,2)
oo <- my.vector.access(v,a)
oo2 <- my.vector.access(v,a,func=prod,zero=1)