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

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

2017-02-27 接線・接円

ryamada2017-02-27

[][][]サインカーブの接線・接円 16:48

  • 接円の中心は、サインカーブの峰のときには鉛直上・下のy=0
  • サインカーブがy=0をよぎるところでは、接円中心は無限遠であり、それは、サインカーブがy=0をよぎる接線に垂直な線方向
  • 途中はこの(pi/2,0) から無限遠を示す直線に向かう漸近曲線

f:id:ryamada:20170226164749j:image

x <- seq(from=0,to=10,length=1000)
y <- sin(x)
t <- seq(from=-pi,to=pi,length=100)
plot(x,y,type="l",xlim=c(-10,10),ylim=c(-10,10))
for(i in 1:length(x)){
	a <- cos(x[i])
	X <- x[i] + t
	Y <- y[i] + a*t
	points(X,Y,type="l",col=2)
	b <- -sin(x[i])
	v <- c(-a,1)
	v <- v/sqrt(sum(v^2))
	ctr <- c(x[i],y[i]) + 1/b * v

	X2 <- ctr[1] + cos(t) * 1/b
	Y2 <- ctr[2] + sin(t) * 1/b
	points(X2,Y2,type="l",col=3)
	points(ctr[1],ctr[2],pch=20,size=1,col=4)
}
トラックバック - http://d.hatena.ne.jp/ryamada/20170227

2017-02-26 京都大学学部入試数学問題2017

[][][]京都大学学部入試数学問題2017で学ぶR 08:29

  • 昨日は国立大学(前期)入学試験の初日。京都大学でも国語と数学の試験がありました。
  • こちらに数学の試験問題(と解答)があります
  • 解説はお任せするとして、Rの練習課題として扱ってみましょう
  • (1)複素数と複素数での変換の話
    • 題意の通りにプロットしてみます

f:id:ryamada:20170226082629j:image

f:id:ryamada:20170226082630j:image

# [1]
## (1)
R <- 2
thetas <- seq(from=0,to=2*pi,length=101)
#thetas <- thetas[-1]

ws <- R * cos(thetas) + R * 1i * sin(thetas)
xy <- ws + 1/ws
xlim <- ylim <- range(c(Re(xy),Im(xy)))

plot(xy,type="l",xlim=xlim,ylim=ylim)

Rs <- seq(from=1,to=2,length=21)
Rs <- Rs[-1]
for(i in 1:length(Rs)){
	R <- Rs[i]
	ws <- R * cos(thetas) + R * 1i * sin(thetas)
	xy <- ws + 1/ws
	xlim <- ylim <- range(c(Re(xy),Im(xy)))
	points(xy,type="l",xlim=xlim,ylim=ylim,col=i+1)

}

## (2)

alpha <- pi/3
Rs <- seq(from=0,to=50,length=1000)
Rs <- Rs[-1]
ws <- Rs * cos(alpha) + Rs * 1i * sin(alpha)
xy <- ws + 1/ws

xlim <- ylim <- range(c(Re(xy),Im(xy)))
xlim[1] <- xlim[1]-5
xlim[2] <- xlim[2]+5
ylim[1] <- ylim[1]-5
ylim[2] <- ylim[2]+5

plot(xy,type="l",xlim = xlim,ylim=ylim)

alphas <- seq(from=0,to=pi/2,length=22)
alphas <- alphas[-1]
alphas <- alphas[-length(alphas)]

for(i in 1:length(alphas)){
	alpha <- alphas[i]
	ws <- Rs * cos(alpha) + Rs * 1i * sin(alpha)
	xy <- ws + 1/ws
	points(xy,type="l",col=i+1)
}
  • (2)3次元の図形。ベクトル。座標の自由度など
    • 平行になること、斜めになること、きれいな正八面体になること

f:id:ryamada:20170226082631p:image

f:id:ryamada:20170226082633p:image

f:id:ryamada:20170226082632p:image

# [2]
# O: X[1,]
# A: X[2,]
# B: X[3,]
# C: X[4,]
# D: r[1] * X[1,] + (1-r[1]) * X[2,]
# E: r[4] * X[2,] + (1-r[4]) * X[3,]
# F: r[6] * X[3,] + (1-r[6]) * X[4,]
# G: r[3] * X[4,] + (1-r[3]) * X[1,]
# H: r[2] * X[1,] + (1-r[2]) * X[3,]
# I: r[5] * X[2,] + (1-r[5]) * X[4,]

X <- matrix(rnorm(3*4),ncol=3)
X <- X/sqrt(apply(X^2,1,sum))
X. <- rbind(X,c(1,1,1))
X. <- rbind(X.,c(-1,-1,-1))
library(rgl)
plot3d(X.,)
spheres3d(X,radius=0.1)
r <- runif(6)
cnt <- 1
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r[cnt] + X[j,] * (1-r[cnt])
		spheres3d(tmp,color=2,radius=0.1)
		cnt <- cnt+1
	}
}

## (1)
r2 <- r
r2[1] <- r2[3]
r2[4] <- 1-r2[6]

cnt <- 1
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
		spheres3d(tmp,color=3,radius=0.1)
		cnt <- cnt+1
	}
}

O <- X[1,]
A <- X[2,]
B <- X[3,]
C <- X[4,]
D <- r2[1] * X[1,] + (1-r2[1]) * X[2,]
E <- r2[4] * X[2,] + (1-r2[4]) * X[3,]
F <- r2[6] * X[3,] + (1-r2[6]) * X[4,]
G <- r2[3] * X[1,] + (1-r2[3]) * X[4,]
H <- r2[2] * X[1,] + (1-r2[2]) * X[3,]
I <- r2[5] * X[2,] + (1-r2[5]) * X[4,]
plot3d(X.,)
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
	}
}
spheres3d(X,radius=0.1)
spheres3d(D,radius=0.1,col=4)
spheres3d(G,radius=0.1,col=4)
spheres3d(E,radius=0.1,col=5)
spheres3d(F,radius=0.1,col=5)

DG <- G-D
EF <- F-E
DG/EF


## (2)

r2 <- rep(1/2,6)

cnt <- 1
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
		spheres3d(tmp,color=3,radius=0.1)
		cnt <- cnt+1
	}
}

O <- X[1,]
A <- X[2,]
B <- X[3,]
C <- X[4,]
D <- r2[1] * X[1,] + (1-r2[1]) * X[2,]
E <- r2[4] * X[2,] + (1-r2[4]) * X[3,]
F <- r2[6] * X[3,] + (1-r2[6]) * X[4,]
G <- r2[3] * X[1,] + (1-r2[3]) * X[4,]
H <- r2[2] * X[1,] + (1-r2[2]) * X[3,]
I <- r2[5] * X[2,] + (1-r2[5]) * X[4,]
plot3d(X.,)
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
	}
}
spheres3d(X,radius=0.1)
spheres3d(D,radius=0.1,col=4)
spheres3d(G,radius=0.1,col=4)
spheres3d(E,radius=0.1,col=5)
spheres3d(F,radius=0.1,col=5)
spheres3d(H,radius=0.1,col=4)
spheres3d(I,radius=0.1,col=4)

simplex <- function(n) {
	qr.Q(qr(matrix(1,nrow=n)),complete=T)[,-1]
}
X <- simplex(4)
r2 <- rep(1/2,6)

cnt <- 1
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
		spheres3d(tmp,color=3,radius=0.1)
		cnt <- cnt+1
	}
}

O <- X[1,]
A <- X[2,]
B <- X[3,]
C <- X[4,]
D <- r2[1] * X[1,] + (1-r2[1]) * X[2,]
E <- r2[4] * X[2,] + (1-r2[4]) * X[3,]
F <- r2[6] * X[3,] + (1-r2[6]) * X[4,]
G <- r2[3] * X[1,] + (1-r2[3]) * X[4,]
H <- r2[2] * X[1,] + (1-r2[2]) * X[3,]
I <- r2[5] * X[2,] + (1-r2[5]) * X[4,]
plot3d(X.,)
for(i in 1:3){
	for(j in (i+1):4){
		segments3d(X[c(i,j),])
		tmp <- X[i,] * r2[cnt] + X[j,] * (1-r2[cnt])
	}
}
spheres3d(X,radius=0.1)
spheres3d(D,radius=0.1,col=4)
spheres3d(G,radius=0.1,col=4)
spheres3d(E,radius=0.1,col=5)
spheres3d(F,radius=0.1,col=5)
spheres3d(H,radius=0.1,col=4)
spheres3d(I,radius=0.1,col=4)
  • (3)整数、三角関数、チョッキリなこと
    • ひたすら調べる

f:id:ryamada:20170226084545p:image

# [3]

p <- q <- 1:20
pq <- as.matrix(expand.grid(p,q))
P <- pq[,1]
Q <- pq[,2]
alpha <- atan(1/P)
beta <- atan(1/Q)
tans <- tan(alpha+2*beta)

plot3d(P,Q,tans)
spheres3d(P,Q,2,col=2,radius=0.1)

pq[which(abs(tans-2)<10^(-12)),]
  • (4)円、三角形、三角関数
    • 外接円を原点中心単位円とし、Aを(1,0)に固定して、BをAから外接円周上、半周させる。それに合わせてCを決める
    • 円の面積を用いて内接円半径を求め、Aから角Aの二等分線上にしかるべき距離の点として内接円中心を描く

f:id:ryamada:20170226092108j:image

f:id:ryamada:20170226092107j:image

# [4]

## 外接円は原点中心の単位円
## A = (1,0)

A <- c(1,0)
thetas <- seq(from=0,to=pi,length=31)
thetas <- thetas[-1]
Bs <- cbind(cos(thetas),sin(thetas))

Cs <- matrix(0,length(thetas),2)
Rs <- rep(0,length(thetas))
InCtr <- matrix(0,length(thetas),2)
for(i in 1:length(thetas)){
	v.B <- Bs[i,]-A
	v.B <- v.B/sqrt(sum(v.B^2)) # 単位ベクトル
	tmp.phi <- Arg(v.B[1]+1i*v.B[2])
	tmp.phi2 <- tmp.phi + pi/3
	r <- -2*cos(tmp.phi2)
	Cs[i,] <- c(r * cos(tmp.phi2) + 1,r*sin(tmp.phi2))
	
	# 内接円半径
	## 三角形面積
	lenAB <- sqrt(sum((Bs[i,]-A)^2))
	lenAC <- sqrt(sum((Cs[i,]-A)^2))
	lenBC <- sqrt(sum((Bs[i,]-Cs[i,])^2))
	area2 <- lenAB * lenAC * sin(pi/3)
	Rs[i] <- area2/(lenAB+lenAC+lenBC)
	tmp.phi3 <- tmp.phi + pi/6
	InCtr[i,] <- A + (Rs[i]/sin(pi/6)) * c(cos(tmp.phi3),sin(tmp.phi3))
	
}

t <- seq(from=0,to=2*pi,length=10)

plot(cos(t),sin(t),pch=20,cex=0.1)

for(i in 1:length(thetas)){
	segments(A[1],A[2],Bs[i,1],Bs[i,2],col=i)
	segments(Bs[i,1],Bs[i,2],Cs[i,1],Cs[i,2],col=i)
	segments(A[1],A[2],Cs[i,1],Cs[i,2],col=i)
}

points(InCtr,col=2,pch=20)

plot(thetas,Rs)
  • (5)微分積分
    • 特定の範囲にひかれた曲線とそれをよぎる直線とに挟まれる領域の面積が直線の傾きによって変化する。面積を最小にする直線は?

f:id:ryamada:20170226093436j:image

f:id:ryamada:20170226093438j:image

# 5
x <- seq(from=0,to=sqrt(2),length=10000)
a <- 0.3
y1 <- x * exp(-x)
y2 <- a * x

plot(x,y1,type="l")
abline(0,a)
abline(v=sqrt(2))


diff.y <- y1-y2
s <- sum(abs(diff.y) * (x[2]-x[1]))

# aを変えてみる
as <- seq(from=0,to=5,length=100)
Ss <- rep(0,length(as))
for(i in 1:length(as)){
	a <- as[i]
	y2 <- a * x
	diff.y <- y1-y2
	Ss[i] <- sum(abs(diff.y) * (x[2]-x[1]))
}

plot(as,Ss,type="l")
abline(h=min(Ss),col=2)


min.a <- as[which(Ss==min(Ss))]
a <- min.a
y1 <- x * exp(-x)
y2 <- a * x

plot(x,y1,type="l")
abline(0,a)
abline(v=sqrt(2))
  • (6)1,2,3,4,5でできたn桁の自然数が3で割り切れる確率
    • ランダム発生させてみる

f:id:ryamada:20170226095604j:image

f:id:ryamada:20170226095605j:image

# まずはべたにやってみる
# [6]

rep.num <- 10^6

max.n <- 10

r <- matrix(sample(1:5,max.n*rep.num,replace=TRUE),ncol=max.n)

frac.mod3 <- rep(0,max.n)
frac.mod3[1] <- length(which(r[,1] %% 3==0))/rep.num

frac.mod3_2 <- frac.mod3
for(i in 2:max.n){
	r10 <- r[,1:i] %*% 10^(0:(i-1))
	tmp <- r10 %% 3
	frac.mod3[i] <- length(which(tmp==0))/rep.num
	tmp.sum <- apply(r[,1:i],1,sum) # これでも同じ
	frac.mod3_2[i] <- length(which(tmp.sum%%3==0))/rep.num 
}


plot(1:max.n,frac.mod3)

plot(frac.mod3,frac.mod3_2)
frac.mod3-frac.mod3_2
# 各桁の値の和の割り切れるかどうかを調べてもよいことを確認

# 桁の値の和の計算の方が楽なので、それを使って桁数を増やしてやってみる
# 桁数を増やすことで、各桁での試行回数を減らしても、挙動の情報は取れるので、そのようにして見る
max.n <- 100
rep.num <- 10^4
r <- matrix(sample(1:5,max.n*rep.num,replace=TRUE),ncol=max.n)

frac.mod3 <- rep(0,max.n)
frac.mod3[1] <- length(which(r[,1] %% 3==0))/rep.num

frac.mod3_2 <- frac.mod3
for(i in 2:max.n){
	# 大変な計算をはしょる
	#r10 <- r[,1:i] %*% 10^(0:(i-1))
	#tmp <- r10 %% 3
	#frac.mod3[i] <- length(which(tmp==0))/rep.num
	tmp.sum <- apply(r[,1:i],1,sum) # これでも同じ
	frac.mod3_2[i] <- length(which(tmp.sum%%3==0))/rep.num 
}

plot(1:max.n,frac.mod3_2)
トラックバック - http://d.hatena.ne.jp/ryamada/20170226

2017-02-20 OpenCV with python

ryamada2017-02-20

[][]Numpy (python)の一番良かった資料 14:34

f:id:ryamada:20170220085849p:image

  • OpenCVのためにnumpyの使い方調べの参照先を決める
  • こちら

[][]OpenCVをとにかく使ってみる 14:44

git clone https://github.com/pmneila/PyMCubes
cd PyMCubes
sudo pip install Cython
sudo python setup.py install
sudo pip install PyCollada
  • その上で ! X,Y,Zとuとの座標xボクセル値の4次元データを、3次元画像データとして渡せばよい
import numpy as np
import mcubes

# Create a data volume (30 x 30 x 30)
X, Y, Z = np.mgrid[:10, :10, :10]
# たくさんの3次元座標に「値」をつける
u = (X-5)**2 + (Y-5)**2 + (Z-5)**2 - 2**2
np.shape(u)
print X
print Y
print u
# Extract the 0-isosurface
# 指定した値0の等高面をマーチングキューブする
vertices, triangles = mcubes.marching_cubes(u, 0)
vertices_xyz = np.transpose(vertices)
print vertices
print triangles
# Export the result to sphere.dae
mcubes.export_mesh(vertices, triangles, "sphere.dae", "MySphere")
  • できたのはColladaフォーマットファイル。これを使うにはPyColladaを使う。その使い方はこちら
from collada import *
mesh = Collada('sphere.dae')
    • プロットする
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
from scipy import genfromtxt

fig = pyplot.figure()
ax = Axes3D(fig)

ax.plot(vertices_xyz[0,],vertices_xyz[1,],vertices_xyz[2,], "o")
col = ["g","r","b"]
cnt = 0

for i,j,k in triangles:
    thiscol = col[cnt]
    ord = [i,j,k,i]

    if(cnt>2):
        cnt = 0
    for t in range(3):
        ax.plot([vertices[ord[t],0],vertices[ord[t+1],0]],[vertices[ord[t],1],vertices[ord[t+1],1]],[vertices[ord[t],2],vertices[ord[t+1],2]],color = col[cnt])
        cnt += 1
        if(cnt>2):
            cnt = 0

pyplot.show()
    • 三角形プロットする
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
import pylab as pl
import scipy as sp

ax = a3.Axes3D(pl.figure())
ax.plot(vertices_xyz[0,],vertices_xyz[1,],vertices_xyz[2,], "o")
cnt = 0
for i,j,k in triangles:       
    vtx = np.array([vertices[i,],vertices[j,],vertices[k,]])
    tri = a3.art3d.Poly3DCollection([vtx])
    tri.set_color(colors.rgb2hex(sp.rand(3)))
    tri.set_edgecolor('k')
    cnt += 1
    #if cnt < 100000:
    ax.add_collection3d(tri)
pl.show()
    • 4D オブジェクトには未対応
import numpy as np
import mcubes

# Create a data volume (30 x 30 x 30)
X, Y, Z, W = np.mgrid[:10, :10, :10, :10]
# たくさんの3次元座標に「値」をつける
u = (X-5)**2 + (Y-5)**2 + (Z-5)**2 - + (W-5)**2  - 2**2
np.shape(u)

# Extract the 0-isosurface
# 指定した値0の等高面をマーチングキューブする
vertices, triangles = mcubes.marching_cubes(u, 0)
      • エラー
RuntimeError                              Traceback (most recent call last)
<ipython-input-29-bad4ae16da1c> in <module>()
     10 # Extract the 0-isosurface
     11 # 指定した値0の等高面をマーチングキューブする
---> 12 vertices, triangles = mcubes.marching_cubes(u, 0)
     13 
     14 print vertices

mcubes/src/_mcubes.pyx in _mcubes.marching_cubes (mcubes/src/_mcubes.cpp:1590)()

RuntimeError: Only three-dimensional arrays are supported.

||<

  • Colladaフォーマットを含むメッシュデータの可視化アプリ meshlabのインストール(こちら)
sudo add-apt-repository ppa:zarquon42/meshlab
sudo apt-get update
sudo apt-get install meshlab
    • 起動する
meshlab
    • 保存したdaeファイルを読み込む
  • pythonのMAYAVIで3Dプロットするにはこちら
トラックバック - http://d.hatena.ne.jp/ryamada/20170220

2017-02-19 PythonをRとの比較で練習する

[][][][]PythonをRとの比較で練習する 10:55

  • こちらのkindleRコードをなぞってみる
  • R notebookではpythonも使える(pythonが普通に動いている環境で〜Windowsとかだと特別な設定が必要、Linuxだと不要
    • LinuxのRにrglをインストールするときのトラブルメモ→こちら
n <- 10
x <- sample(1:3,n,replace=TRUE)
x
mean(x)
1/n * sum(x)

x
tabulate(x)
w <- tabulate(x)/n
w
v <- sort(unique(x))
v
mw <- sum(v * w)
mw
>|python|
n = 10
import numpy as np

x = np.array(np.random.randint(1,4,n),dtype="f") # replace=TRUE is default
print x
print np.average(x)
print np.sum(x)/n
import collections
tab = collections.Counter(x)
print tab
m_w = 0.
for k,v in tab.items():
  print k,v
  m_w += k*v
print m_w/len(x)
v <- 1:6
p <- rep(1/6,6)
sum(v*p)

v <- 1:6
p <- rep(1/6,6) + rnorm(6) * 0.1
p <- p/sum(p)
plot(v,p,type="h")
p
sum(p)
sum(v*p)
import numpy as np
v = np.arange(1,7,dtype="f")
print v
p = np.repeat([1./6],6)
print p
print sum(v*p)

p += np.random.randn(6)*0.1
p /= np.sum(p)
print p
print np.sum(p)
print np.sum(v*p)
import matplotlib.pyplot as plt
#%pylab inline --no-import-all
plt.bar(np.arange(1,7),p)
#plt.show()
p0 <- 0.3
n <- 10
i <- 0:n
i.inv <- i[(n+1):1]
choose(n,i)
p <- choose(n,i) * p0^i * (1-p0)^i.inv
plot(i,p,type="h")
sum(p)
sum(p*i)

n * p0
import numpy as np
from scipy import special
p0 = 0.3
n = 10
i = np.arange(n+1)
print i
i_inv = i[::-1]
print i_inv
p = np.zeros(n+1)
for k in range(n+1):
  p[k] = special.binom(n,k) * p0**i[k] * (1-p0)**i_inv[k]

print p 
print special.binom(n,i) * p0**i * (1-p0)**i_inv
p = special.binom(n,i) * p0**i * (1-p0)**i_inv
print np.sum(p)
print np.sum(p*i)
p <- seq(from=0,to=1,length=100)
h <- p
plot(p,h,type="l")
import numpy as np
n = 100
p = np.arange(n+1,dtype="f")/n
print p
n <- 3
m <- 5
p <- 0.2
gamma(n+m+2)/(gamma(n+1)*gamma(m+1))*p^n*(1-p)^m

dbeta(p,n+1,m+1)

p <- seq(from=0,to=1,length=100)
d <- dbeta(p,n+1,m+1)
plot(p,d,type="l")
n = 3
m = 5
p = 0.2
from scipy import stats, special
print special.gamma(n+m+2)/(special.gamma(n+1)*special.gamma(m+1)) * p**n * (1-p)**m
print stats.beta.pdf(p,n+1,m+1)
import numpy as np
p = np.arange(101,dtype="f")/100
d = stats.beta.pdf(p,n+1,m+1)
#print d
import matplotlib.pyplot as plt
plt.plot(p,d)
plt.vlines(n/(n+m+0.),0,np.max(d))
plt.show()
lambda <- 2.8
n <- 1000
x <- rpois(n,lambda)
plot(table(x))
mean(x)
from scipy import stats
k = 2.8
n = 1000
x = stats.poisson.rvs(k,size=n)
#print x
import collections
tab = collections.Counter(x)
print tab
ns <- 0:20
p <- lambda^ns/factorial(ns) * exp(-lambda)
plot(ns,p,type="h")

plot(ns,dpois(ns,lambda),type="h")
import numpy as np
from scipy import misc, stats
ns = np.arange(21)
k = 2.8
p = k**ns/misc.factorial(ns) * np.exp(-k)
print p
import matplotlib.pyplot as plt
plt.bar(ns,p)
plt.show()
n <- 4
lambdas <- seq(from=0,to=10,length=100)
L <- lambdas^n/factorial(n) * exp(-lambdas)
plot(lambdas,L,type="l")
plot(lambdas,log(L),type="l")
import numpy as np
from scipy import special
n =4
k = np.arange(101)/100.*10
L = k**n/special.factorial(n) * np.exp(-k)
import matplotlib.pyplot as plt
plt.plot(k,L)
plt.show()
plt.plot(k,np.log(L))
plt.show()

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

2017-02-12 pythonの練習

[][][][]numpy,scipy,sympyを実務に使って練習する 15:53

  • ライブラリ・モジュール、その読み込み
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pyshtools as sh
from sympy.physics.quantum.cg import CG,cg_simp
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import Delaunay
from scipy.special import sph_harm
import math
import collections # 多値返却のため
%matplotlib notebook
    • import hoge as hg は、呼び出しを簡単にしてhg.hageのようにして使うため
    • from hoge import hage はhogeの中からhageを特に読み込んで、hageをそのままの名前で使うため
  • ベクトル・アレイを使って三角形の計算する
    • 重心
# 三角形の重心座標を三角化オブジェクトから作る
def triweightedcenter(ts):
    ret = np.zeros((3,ts.tri.T[0,].size))
    for i in range(ts.tri.T[0,].size):
        id1 = ts.tri[i,0]
        id2 = ts.tri[i,1]
        id3 = ts.tri[i,2]

        v1 = np.array([ts.X[0,id1],ts.X[1,id1],ts.X[2,id1]])
        v2 = np.array([ts.X[0,id2],ts.X[1,id2],ts.X[2,id2]])
        v3 = np.array([ts.X[0,id3],ts.X[1,id3],ts.X[2,id3]])
        ret[0,i], ret[1,i], ret[2,i] = (v1 + v2 + v3)/3
        
    return ret
    • 頂点座標を持つ2次元アレイと頂点IDの三つ組を渡して、複数の三角形の重心を計算する(頂点IDトリオ、頂点座標2次元アレイのハンドリング関数は後述)
def triweightedcenter2(tri,X):
    ret = np.zeros((3,tri.T[0,].size))
    for i in range(tri.T[0,].size):
        id1 = tri[i,0]
        id2 = tri[i,1]
        id3 = tri[i,2]

        v1 = np.array([X[0,id1],X[1,id1],X[2,id1]])
        v2 = np.array([X[0,id2],X[1,id2],X[2,id2]])
        v3 = np.array([X[0,id3],X[1,id3],X[2,id3]])
        ret[0,i], ret[1,i], ret[2,i] = (v1 + v2 + v3)/3
        
    return ret
    • 面積
# 三角形の面積を計算する
# 面積は |u||v|sin(theta)/2
# v1,v2,v3はnp.arrayなので、ベクトル演算(要素ワイズな加減乗除など)ができるし、内積に相当するdot計算もできる
def triarea(v1,v2,v3):
    e12 = v2-v1
    e13 = v3-v1
    ip = np.dot(e12,e13)
    r1 = np.sqrt(np.sum(np.power(e12,2)))
    r2 = np.sqrt(np.sum(np.power(e13,2))) 
    costheta = ip/(r1*r2)
    sintheta = np.sqrt(1-np.power(costheta,2))
    area = 0.5 * r1*r2*sintheta
    return area
    • 頂点座標を持つ2次元アレイと頂点IDの三つ組を渡して、複数の三角形の面積を計算する
# 三角形の面積を頂点ID三つ組リストと頂点座標格納行列とから作る
# 頂点IDは、三角形の個数だけ、長さ3の整数型リストを束ねた2次元アレイ
# 頂点座標は、頂点の個数の長さのベクトルを3つ束ねたもの
def triarea_mult(tri,X):
    ret = np.zeros(tri.T[0,].size)
    for i in range(tri.T[0,].size):
        id1 = tri[i,0]
        id2 = tri[i,1]
        id3 = tri[i,2]
        if (id1-id2)*(id2-id3)*(id3-id1) == 0:
            ret[i] = 0
        else:
            v1 = np.array([X[0,id1],X[1,id1],X[2,id1]])
            v2 = np.array([X[0,id2],X[1,id2],X[2,id2]])
            v3 = np.array([X[0,id3],X[1,id3],X[2,id3]])
            ret[i] = triarea(v1,v2,v3)
        
    return ret    

    • 三角形メッシュを描図する
# 三角化プロット
def spheretriplot(ts):  
# 3Dplot
    fig = plt.figure() # 描く領域を確保して
    ax = Axes3D(fig)
# The triangles in parameter space determine which x, y, z points are
# connected by an edge
# 頂点と三角形エッジを描く
# 色具合のオプションを使って少しきれいに
    ax.plot_trisurf(ts.X[0,],ts.X[1,],ts.X[2,], triangles=ts.tri, cmap=plt.cm.Spectral)
# 表示させる
    plt.show()
  • 球面を覆う三角形メッシュを作る
# 経線、緯線数を与えて三角グリッドを作る
# nlatは緯線の数、nlat = 3は、北極、赤道、南極の3つ
# nlonは経線の数、nlon = 2は、グリニッヂと日付変更線の2つ
# 三角形は、2本の緯線の間に蛇腹形式で、北極・南極は、隣接緯線と多角錐をなす
# phi,thetaは極座標系 https://en.wikipedia.org/wiki/Polar_coordinate_system
def spheregrid(nlat=10,nlon=10):
    phi = 2*math.pi/nlon * np.array([range(nlon)])
    theta = math.pi/(nlat-1) * np.array([range(nlat)])
    
    n = 2 + nlon * (nlat-2) # 頂点数
    
    ret = np.zeros((2,n)) # phi,thetaの格納庫
    ret[0,n-1] = math.pi
    for i in range(nlat -2):
        for j in range(nlon):
            tmp = i * nlon + j + 1
            ret[0,tmp] = theta[0,i+1]
            ret[1,tmp] = phi[0,j]
    
    n_tri = nlon * 2 + nlon*2*(nlat-3) # 三角形の数
    tri = np.zeros((n_tri,3),dtype=np.int) # 三角形の3頂点idの格納庫、idは整数型
    cnt = 0
    for i in range(nlon-1):
        tri[cnt,0] = 0
        tri[cnt,1] = i+1
        tri[cnt,2] = i+2
        cnt +=1
    tri[cnt,0] = 0
    tri[cnt,1] = nlon
    tri[cnt,2] = 1
    cnt +=1
    for i in range(nlat-3):
        for j in range(nlon-1):
            tri[cnt,0] = i*nlon + j + 1
            tri[cnt,1] = i*nlon + j + 2
            tri[cnt,2] = (i+1)*nlon + j + 2
            cnt += 1
            tri[cnt,0] = i*nlon + j + 1
            tri[cnt,1] = (i+1)*nlon + j + 1
            tri[cnt,2] = (i+1)*nlon + j + 2
            cnt += 1
        tri[cnt,0] = i*nlon + nlon
        tri[cnt,1] = i*nlon + 1
        tri[cnt,2] = (i+1)*nlon + 1
        cnt += 1
        tri[cnt,0] = i*nlon + nlon
        tri[cnt,1] = (i+1)*nlon + nlon
        tri[cnt,2] = (i+1)*nlon + 1
        cnt += 1
    for i in range(nlon-1):
        tri[cnt,0] = (nlat-3)*nlon + i + 1
        tri[cnt,1] = (nlat-3)*nlon + i + 2
        tri[cnt,2] = n-1
        cnt +=1
    tri[cnt,0] = (nlat-3)*nlon + nlon-1 + 1
    tri[cnt,1] = (nlat-3)*nlon + 1
    tri[cnt,2] = n-1
    
    X = np.ndarray((3,n)) # 頂点のデカルト座標(原点中心単位球の場合)
    for i in range(ret[0,].size):
        X[0,i], X[1,i], X[2,i] = np.sin(ret[0,i])*np.cos(ret[1,i]), np.sin(ret[0,i])*np.sin(ret[1,i]), np.cos(ret[0,i])
    # ここまでで頂点座標情報と、三角形頂点ID情報が作られてる
  # 以下はこの球面三角メッシュの属性の計算
    # 三角形の面積
    areas = triarea_mult(tri,X)
    # 三角形の重心座標
    weightedcenter = triweightedcenter2(tri,X)
    #print weightedcenter
    # 頂点と三角形重心の極座標(後述)
    #Xpolar = decart2polar3d_ndarray(ts.X)
    WCpolar = decart2polar3d_ndarray(weightedcenter)
    
    result = collections.namedtuple('result', 'X, vangles, fangles,tri,area,wcenter,')
    return result(X=X, vangles=ret, fangles=WCpolar, tri=tri, area = areas, wcenter=weightedcenter)
  • 3D 極座標系
# 3Dデカルト座標を極座標に変える(r==0は考えない)
def decart2polar3d(x,y,z):
    r = math.sqrt(x**2+y**2+z**2)
    rxy = x**2+y**2
    if rxy==0:
        phi = 0
        theta = 0
        if z < 0:
            theta = math.pi
    else:
        if x==0:
            phi = math.pi/2
            if y < 0:
                phi = -math.pi/2
        else:
            phi = np.arctan(y/x)
        theta = np.arccos(z/r)
    return [phi,theta,r]
    
decart2polar3d(0,0,-1)

def decart2polar3d_ndarray(X):
    ret = np.zeros_like(X)
    for i in range(X[0,].size):
        tmp= decart2polar3d(X[0,i],X[1,i],X[2,i])
        ret[0,i],ret[1,i],ret[2,i] = tmp
    return ret
  • 球面場を実数球面調和関数の和として考えることにする
    • 球面調和関数係数をpyshtoolsのSHCoeffsクラスのcoeffsメンバーに登録する形で定める
    • ランダムに与えることにする
    • degree l までの係数は (l+1)^2個となるが、それを(l+1)^2 * 2個の要素を格納できる3次元アレイに格納するのがSHCoeffsクラスのやりかた
# 球面になめらかな場を与える
# 球面調和関数の係数のセットとして与える
# x.coeffsは[0,,]にm>=0[1,,]にm<0が格納される
# [*,0,k]はabs(m)=kの値を格納する
def randSHcoef(lmax):
    num_coef = (lmax+1)**2
    value_re = np.random.randn(num_coef) + 0j
#print value_re
    x = sh.SHCoeffs.from_zeros(lmax,'complex')

    for l in range(0, lmax+1):
        for m in range(-l, l+1):
          x.set_coeffs(value_re[l*(l+1)+m], l, m)
    return x
  • 球面場の値を極座標に与える
# scipy.special.sph_harmを使ってSHcoeffから、球面座標の値を計算
def coeff2valreal(cf,phi,theta):
    ret = 0
    #lmax = cf
    for l in range(0, cf[0,0,].size):
        for m in range(-l, l+1):
            if m >= 0:
                tmp = cf[0,l,m]
                ret += sph_harm(m,l,phi,theta).real
            else:
                tmp = cf[1,l,np.abs(m)]
                ret += sph_harm(m,l,phi,theta).imag
    return ret

def coeff2valcomplex(cf,phi,theta):
    ret = 0+1j
    #lmax = cf
    for l in range(0, cf[0,0,].size):
        for m in range(-l, l+1):
            if m >= 0:
                tmp = cf[0,l,m]
            else:
                tmp = cf[1,l,np.abs(m)]
            ret += sph_harm(m,l,phi,theta)
    return ret

# 三角の頂点・重心座標の場の値を球面調和関数係数numpy.ndarrayから作る
def SHval(Xpolar,coef):
    valreal = np.zeros(Xpolar[0,].size)
    valcomplex = np.zeros(Xpolar[0,].size,'complex')
    for i in range(Xpolar[0,].size):
    #valreal += [ coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])]
    #valcomplex += [coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])]
        valreal[i] = coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])
        valcomplex[i] = coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])
        
    result = collections.namedtuple('result', 'real, clx')
    return result(real = valreal,clx=valcomplex)
  • 球面調和関数分解、その離散版
    • 個々の球面調和関数の分解係数は、球面全体にわたって、場の関数と球面超関数との積分を取ることである
    • 離散版は、小領域ごとに場関数の値と球面調和関数の値との積を取り、それを領域全体にわたって、小領域面積で重み付けした積算を取ること
    • その処理を個々の球面調和関数についておこない、それを指定の球面調和関数セットについて行う
# 球面調和関数分解する
# ある球面調和関数の係数は:
# 全三角形にて総和を取る(三角形重心の値 x 三角形重心のある球面調和関数の値 x 三角形の面積)
# 実関数版球面調和関数分解の場合には、m>=0は実部を、m<0は虚部をとる。それは位相をpi/2ずらすことに相当
# 複素環数版分解の場合はそのまま。そうすると m,-mとの係数は共役複素数になる
def SHdecomp(fvalue,area,fangle,l,m):
    ret_re = 0
    ret_clx = 0 + 0j
        
    for i in range(area.size):
        if m < 0:
            ret_re += fvalue[i] * sph_harm(m,l,fangle[0,i],fangle[1,i]).imag * area[i]
        else:
            ret_re += fvalue[i] * sph_harm(m,l,fangle[0,i],fangle[1,i]).real * area[i]
        ret_clx += fvalue[i] * sph_harm(m,l,fangle[0,i],fangle[1,i]) * area[i]
    return [ret_re,ret_clx]

# lmaxまでのすべての分解係数を3次元アレイに格納する
# [2,l+1,l+2]アレイで[0,,]はm>=0,[1,,]はm<0
def SHdecomp_lmax(fvalue,area,fangle,lmax):
    xre = sh.SHCoeffs.from_zeros(lmax)
    xclx = sh.SHCoeffs.from_zeros(lmax,'complex')
    for l in range(0, lmax+1):
        for m in range(-l, l+1):
            tmp = SHdecomp(fvalue,area,fangle,l,m)
            xre.set_coeffs(tmp[0], l, m)
            xclx.set_coeffs(tmp[1], l, m)
    result = collections.namedtuple('result', 'real, complex')
    return result(real = xre,complex=xclx)
  • やってみる
ts = spheregrid(nlat=10,nlon=10)
spheretriplot(ts)
# 球面場モデルを作る
f = randSHcoef(50)
# 場の値を計算
fvalV = SHval(ts.vangles,f.coeffs).real
fvalF = SHval(ts.fangles,f.coeffs).real
l = 2
m = 1
lmax = 10
fakefvalF = np.random.randn(ts.area.size)
print SHdecomp(fakefvalF,ts.area,ts.fangles,l,m)
dcreal = SHdecomp_lmax(fakefvalF,ts.area,ts.fangles,lmax).real
dcclx = SHdecomp_lmax(fakefvalF,ts.area,ts.fangles,lmax).complex
print dcreal.coeffs
print dcclx.coeffs
alpha = math.pi/6
beta = math.pi/4
gamma = math.pi/3
  • これのよいところは、球を回転させた後の球面調和関数係数を、x.rotate(a,b,c)とするだけでpyshtoolsの関数で計算できること
# 係数はそのまま回転できる
dcreal_rot = dcreal.rotate(alpha,beta,gamma)
dcclx_rot = dcclx.rotate(alpha,beta,gamma)
print dcreal_rot.coeffs
print dcclx_rot.coeffs
  • 回転によって、分解係数の二乗ノルムは変わらない。また、lごとの係数二乗ノルムも変わらない
    • それを確かめる
# 実SH分解では、回転前後でlごとの係数ベクトル二乗ノルムが変わらないこと
print "REAL"
print "all"
print np.sum(np.power(dcreal.coeffs,2))
print np.sum(np.power(dcreal_rot.coeffs,2))
print "--"
for i in range(lmax+1):
    print i
    print np.sum(np.power(dcreal.coeffs[0,i,],2)+np.power(dcreal.coeffs[1,i,],2))
    print np.sum(np.power(dcreal_rot.coeffs[0,i,],2)+np.power(dcreal_rot.coeffs[1,i,],2))
print "--"
# 複素SH分解では二乗ノルムが(a+bj)*(a-bj)であることに注意
# 実 a+bj -> m>0 : a. m<0: b -> a^2 + b^2
# 複 (m,-m) : (a+bj, a-bj) -> (a+bj)(a-bj) + (a-bj)(a+bj) = (a^2 + b^2)*2
# (a+bj)(a-bj)は、各係数の長さの二乗でもあるが、clmとcl(-m)の積でもあるので、この値の意味合いは少し変わってくる
# その先にrotation invariantsがあるような気がする
print "COMPLEX"
print "all"
print np.sum(np.power(abs(dcclx.coeffs),2))
print np.sum(np.power(abs(dcclx_rot.coeffs),2))
print "--"
for i in range(lmax+1):
    print i
    print np.sum(np.power(abs(dcclx.coeffs[0,i,]),2)+np.power(abs(dcclx.coeffs[1,i,]),2))
    print np.sum(np.power(abs(dcclx_rot.coeffs[0,i,]),2)+np.power(abs(dcclx_rot.coeffs[1,i,]),2))
  • sympyによる球面調和関数の回転変換
    • 球面調和関数で表されるものには物理・量子物理の電子の状態などがある。それもあってscipyに球面調和関数のハンドリング関数があるのだが、from sympy.physics.quantum にも、そんなのがある
    • sympyは代数処理(変数を持った式を扱う仕組み。いちいち数値計算しないで、式変形をする仕組み。必要ならば、式変形の最後に式を評価する)
    • 分解係数の回転にあたってClebsch-Goradn係数というのがある。これは係数の回転変換に使うWignerDの相方のようなもので、階乗のおばけのようなもので面倒くさいし、うまくすれば割り算などが式変形でなくせる
    • それを使った関数を作って見る
# Clebsch-Gordan coefficientsのsympyを使う
def vl1l2jk(l1,l2,j,k,coef):
    ret = 0
    for m in range(-l1,l1+1):
        if np.abs(k-m) <= l2:
            if m < 0:
                tmp1 = coef[1,l1,-m]
            else:
                tmp1 = coef[0,l1,m]
            if k-m < 0:
                tmp2 = coef[1,l2,-k+m]
            else:
                tmp2 = coef[0,l2,k-m]
            ret += CG(l1,m,l2,k-m,j,k) * tmp1 * tmp2
        else:
            ret += 0
    return ret

def vl1l2vl3l4j(l1,l2,l3,l4,j,coef):
    ret =0
    for k in range(-j,j+1):
        ret += (-1)**(j-k)*vl1l2jk(l1,l2,j,k,coef)*vl1l2jk(l3,l4,j,-k,coef)
    ret /= np.sqrt(2*j+1)
    return ret

def vl1l2jv(l1,l2,j,coef):
    ret = 0
    for k in range(-j,j+1):
        if (-k) <0:
            tmp = coef[1,l2,k]
        else:
            tmp = coef[0,l2,-k]
        ret += (-1)**(j-k)*vl1l2jk(l1,l2,j,k,coef)*tmp
    ret /= np.sqrt(2*j+1)
    return ret
    • 使ってみる
l1 = 1
l2 = 1
j = 0
k = 0
# sympy の代数形式
print vl1l2jk(l1,l2,j,k,dcclx.coeffs)
print vl1l2jk(l1,l2,j,k,dcclx_rot.coeffs)
# その値算出
print vl1l2jk(l1,l2,j,k,dcclx.coeffs).doit().evalf()
print vl1l2jk(l1,l2,j,k,dcclx_rot.coeffs).doit().evalf()

# 実数係数は不変ではない
print vl1l2jk(l1,l2,j,k,dcreal.coeffs)
print vl1l2jk(l1,l2,j,k,dcreal_rot.coeffs)
print vl1l2jk(l1,l2,j,k,dcreal.coeffs).doit().evalf()
print vl1l2jk(l1,l2,j,k,dcreal_rot.coeffs).doit().evalf()
トラックバック - http://d.hatena.ne.jp/ryamada/20170212