ryamadaのコンピュータ・数学メモ RSSフィード

統計学・遺伝学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
数学用語集(TOMAC)

2018-05-16 my circos

ryamada2018-05-16

[][]my circos 02:04

  • 要素を円周上に並べ、関連のある要素を円の内側に弧を描くデータ視覚化手法がある。circos
  • きれいだけれど、簡易に描きたい
  • 特に、データ構造のみの情報で描きたい
  • このようなデータ構造は色々考えられるが、自分の用途は次のようなもの
    • 木グラフがある。木の周りをぐるりと回ると、それは周回グラフになる。それを円周に配置する
    • 木の枝を辿って根から遠い方に進み、先端で折り返して戻ってきて、枝の生え際に帰ってくることになるが、この枝の生え際を複数回通過することを、「同一視」として弧で結ぶことにする。これにより、木の周回というデータ構造を表す用途に使える
  • 今、周回グラフがあるとする。それは円周上に等間隔で配置できる
  • 周回グラフ上のところどころに、対になるものがあったときに、その間に弧を描きたい
  • たとえば (0,4,1,2,10,3,10,1,4,5,6,7,8,4,9) というような要素列で1周とすると、途中に1が2回、10が2回、4が3回登場する。
  • 今、単位円周上の2点が、角 ¥alpha ¥le ¥betaで指定されているとき、その2点間に弧を引きたい
  • その弧は単位円周と直交するものとすれば、以下のようにその弧は指定できる
    • ¥theta = ¥beta-¥alphaとする。¥frac{¥alpha+¥beta}{2}=¥phiと置く(これは2点の中間方向)
    • 弧を含む円の半径は¥tan{¥frac{¥theta}{2}}
    • 弧を含む円の中心までの距離はr = ¥frac{1}{¥cos{¥frac{¥theta}{2}}}
    • 弧を含む円の中心の座標は(r ¥cos{¥phi}, r¥sin{¥phi})
    • 弧は、角¥frac{¥pi}{2} + ¥betaから¥frac{3¥pi}{2} + ¥alphaまで
  • Rで関数にする
# Yamada-cirq

my.cirq <- function(x,N=1000){
	n <- length(x)
	theta <- (0:(n-1))/n * 2*pi
	X <- cbind(cos(theta),sin(theta))
	d <- as.matrix(dist(x))
	diag(d) <- 1
	pairs <- which(d==0,arr.ind=TRUE)
	pairs <- unique(t(apply(pairs,1,sort)))

	alphas <- theta[pairs[,1]]
	betas <- theta[pairs[,2]]
	
	rs <- 1/cos((betas-alphas)/2)
	
	ctrs <- rs * cbind(cos((alphas+betas)/2),sin((alphas+betas)/2))
	
	Rs <- tan((betas-alphas)/2)
	
	phis <- cbind(betas+pi/2, alphas+3*pi/2)
	Theta <- (0:(N-1))/N * 2 * pi
	plot(cos(Theta),sin(Theta),type="l")
	points(X,pch=20)
	for(i in 1:length(pairs[,1])){
		Phi <- seq(from = phis[i,1],to=phis[i,2],length=N)
		arc <- Rs[i] * rbind(cos(Phi),sin(Phi)) + ctrs[i,]
		points(t(arc),type="l")
	}
	
}

x <- c(0,4,1,2,10,3,10,1,4,5,6,7,8,4,9)

my.cirq(x)

f:id:ryamada:20180517020102j:image

トラックバック - http://d.hatena.ne.jp/ryamada/20180516

2018-05-15 部分集合要素のソート

[][]部分集合要素のソート 11:14

# ある集合に順序があり
# その部分集合を全体集合に定まった順序でソートしたい

my.sort.subset <- function(subset,setorder){
	loc <- rep(0,length(subset))
	for(i in 1:length(loc)){
		loc[i] <- which(setorder==subset[i])
	}
	ord <- order(loc)
	return(subset[ord])
}

a <- c(7,2,5)
b <- sample(1:8)
a
b

my.sort.subset(a,b)
トラックバック - http://d.hatena.ne.jp/ryamada/20180515

2018-05-12 重複なし行のみの抽出

[]重複なし行のみの抽出 22:28

  • データレコードが行列状になっているときに、重複行を見つけて、ダブっている分を除いて、1行分にして出すのがunique()
  • 今、ダブっているレコードは1行残らず、破棄したい
  • unique()関数は、上からスキャンして、初出はTRUE、既出分とダブったらFALSEを返す関数 duplicate()を基準にしているので、スキャンを上下の2方向からやって、どちらでもTRUEな行は、重複に関わらない行
my.onlyone <- function(x){
	n <- length(x[,1])
	a <- !duplicated(x)
	b <- !duplicated(x[n:1,])[n:1]
	return(x[a & b,])
}
x <- matrix(c(0,0,0,1,0,1,1,1),byrow=TRUE,ncol=2)
my.onlyone(x)
  • これを使って、『図形』がブロックで出来ているときに、内部点・内部面などを取り除いて表面要素だけにすることができるだろう
# ボクセルはその座標最小頂点(x,y,z)座標で登録する

# 各ボクセルは6面を持つ
# 面(正方形)はその座標最小頂点(x,y,z)と、面の方向で区別する
# xy平面、yz平面、zx平面に広がるものを、1,2,3と呼び分ける

# 各頂点は4頂点を持つ

Vox.list <- as.matrix(expand.grid(0:2,0:2,0:2))

my.vox.faces <- function(xyz){
	rbind(c(xyz,1),c(xyz,2),c(xyz,3),c(xyz+c(1,0,0),2),c(xyz+c(0,1,0),3),c(xyz+c(0,0,1),1))
}

my.face.nodes <- function(xyzw){
	xyz <- xyzw[1:3]
	w <- xyzw[4]
	if(w==1){
		ret <- rbind(xyz,xyz+c(1,0,0),xyz+c(0,1,0),xyz+c(1,1,0))
	}else if(w==2){
		ret <- rbind(xyz,xyz+c(0,1,0),xyz+c(0,0,1),xyz+c(0,1,1))
	}else{
		ret <- rbind(xyz,xyz+c(0,0,1),xyz+c(1,0,0),xyz+c(1,0,1))
	}
	return(ret)
}

my.onlyone <- function(x){
	n <- length(x[,1])
	a <- !duplicated(x)
	b <- !duplicated(x[n:1,])[n:1]
	return(x[a & b,])
}

my.vox.surface <- function(Vlist){
	faces <- matrix(0,0,4)
	for(i in 1:length(Vlist[,1])){
		faces <- rbind(faces,my.vox.faces(Vlist[i,]))
	}
	return(my.onlyone(faces))
	
}

faces <- my.vox.surface(Vox.list)

my.surface.nodes <- function(faces){
	vs <- matrix(0,0,3)
	for(i in 1:length(faces[,1])){
		vs <- rbind(vs,my.face.nodes(faces[i,]))
	}
	return(unique(vs))
}

nodes <- my.surface.nodes(faces)

library(rgl)
plot3d(nodes)


トラックバック - http://d.hatena.ne.jp/ryamada/20180512

2018-04-29 球面調和関数同士の関係

[]球面調和関数同士の関係 17:59

  • 2つの球面調和関数を掛け合わせたり、¥cos{¥theta}を掛けたりするとどうなるかに関する情報の載っている資料:こちら
トラックバック - http://d.hatena.ne.jp/ryamada/20180429

2018-04-28 pyperを使ってみる

[][]pyperを使ってみる 04:16

  • pythonからRを動かしてその結果をpython内に取り戻すことができる
  • 参考
  • pyperパッケージをインストールして使えるようにしたうえで
  • numpyオブジェクトを使ってRのオブジェクトを取り出すようなRの動かし方を指定する
from pyper import *
r = R(use_numpy=True)
  • Rの1文コマンドを実行して、その結果を取り出す
outputs = r("a <- 3; print(a + 5)")
outputs
  • Rのコンソールに標準出力される内容を吐く
'try({a <- 3; print(a + 5)})\n[1] 8\n'
  • pythonにRオブジェクトを取り込む
x_py =r.get('a')
x_py
  • aに格納されていた 3 が取り出される
3
  • 行列も取り出せる
r("x = matrix(rnorm(20),ncol=2)")
x_py =r.get('x')
x_py
array([[ 0.48296985, -1.31237143],
       [ 0.26313735,  0.17361424],
       [-0.62325386, -1.57476065],
       [-2.2620153 ,  1.83864163],
       [ 1.31291217,  1.59120925],
       [-1.72555325, -0.69653969],
       [ 0.94822337,  1.76162395],
       [-0.0518286 , -0.15631062],
       [ 1.58084999,  0.64470825],
       [ 2.0807602 , -0.07323586]])
  • 長いコードを与えるのは、ちょっと面倒くさい模様。改行を許さないので、";"区切りで連続書きする
rscript = "x = matrix(rnorm(210),ncol=2);"
rscript += "x = x/apply(x^2,1,sum);"
rscript += "obj = list();"
rscript += "library(onion);"
rscript += "qtn = 3+2*Hi;"
rscript += "obj[[1]] = qtn"
rscript
Out[8]:
'x = matrix(rnorm(210),ncol=2);x = x/apply(x^2,1,sum);obj = list();library(onion);qtn = 3+2*Hi;obj[[1]] = qtn'
  • 取り出そう
r("x = " +rscript)
x_py =r.get('obj[[1]]')
x_py
array([[3],
       [2],
       [0],
       [0]])
トラックバック - http://d.hatena.ne.jp/ryamada/20180428