球面調和関数をrPython でやる - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-07-19

球面調和関数をやりたいがRではいいものがないのでPythonでやった結果をrPython を使ってRでゴリ押し計算するやり方。R を投げ捨ててPython に移行したかったが、結局R でやっちゃう。
rPython について知りたければ下にスクロール。
　
球面上の離散的な分布をなんらかの関数で級数展開できれば、関数化できる。球面上の分布 $f$ は、球座標表示で $f(\theta, \phi)$ と書ける。これを、ルジャンドル多項式
$P_l(x)=\frac{1}{2^n n!}\frac{d^l}{dx^l}(x^2-1)^l$
とルジャンドル陪関数
$P_l^m(x)=(-1)^m(1-x^2)^{\frac{m}{2}}\frac{d^m}{dx^m}P_l(x)$
を用いて、球面調和関数
$Y_l^m(\theta, \phi)=(-1)^{\frac{m+|m|}{2}}\sqrt{\frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(\cos\theta)exp(im\phi)$
と表すことができる。虚数が入っているので、実数版は
$Y_l^m=\{\begin{matrix}\sqrt{2}(-1)^mRe(Y_l^m)&if\hspace{3}m>0\\Y_l^0&if\hspace{3}m=0\\\sqrt{2}(-1)^{-m}Im(Y_l^{-m})&if\hspace{3}m<0\end{matrix}$
となる。
l はdegree, m はorder とよく書かれている。 $-l\leq m \leq l$ である。
　
Rでは球面調和関数を直接やってくれそうな関数がなさそうなので、組み合わせてやってみる。orthopolynom パッケージはルジャンドル多項式を返す。結果は0 degree からl degree まで出力される。

library(orthopolynom)
l_spharm <- 5 # l; degree
m_spharm <- 0 # m; order
legp <- legendre.polynomials(l_spharm)

[[1]]
1 

[[2]]
x 

[[3]]
-0.5 + 1.5*x^2 

[[4]]
-1.5*x + 2.5*x^3 

[[5]]
0.375 - 3.75*x^2 + 4.375*x^4 

[[6]]
1.875*x - 8.75*x^3 + 7.875*x^5

legendre.polynomials は直接微分してくれる。複数回の微分(m) はその都度loop を回す。

l_spharm <- 5 # l; degree
m_spharm <- 2 # m; order
legp <- legendre.polynomials(l_spharm)[l_spharm + 1]
for(mm in 1:m_spharm){
  legp <- polynomial.derivatives(legp)
}

[[1]]
-52.5*x + 157.5*x^3

さてここで、具体的なルジャンドル倍関数を求めて確かめてみる。
$P_2^2(x)=-3(1-x^2)$ なので、

l_spharm <- 2 # l; degree
m_spharm <- 2 # m; order
legp <- legendre.polynomials(l_spharm)[l_spharm + 1]
for(mm in 1:m_spharm){
  legp <- polynomial.derivatives(legp)
}
P_l_m <- function(x) sum((-1)^m_spharm * (1-x^2)^(m_spharm/2) *polynomial.coefficients(legp)[[1]]* x^(0:(l_spharm-m_spharm)))

x <- seq(-1, 1, length=100)
plot(x, mapply(P_l_m, x))
#lines(x, -3*x*sqrt(1-x^2)) # l <- 2; m <- 1
lines(x, 3*(1-x^2))

確かにあってる。

　
本当なら、球面調和関数同士の積、内積がクロネッカーのデルタ $\delta_{l,l^{'}}$ になることを利用して、もとの球面上の関数 $f(\theta,\phi)$ は
$f(\theta,\phi)=\sum_{l=0}^\infty\sum_{m=-l}^lc_l^mY_l^m(\theta,\phi)$
$c_l^m=\int_{-{\pi}}^{\pi} d\phi\int_0^{\pi} d\theta \sin\theta f(\theta, \phi)Y_l^m(\theta, \phi)$
となる。
が、全然うまくRで再現できない。
上のリンクでわざわざPython コードがあるので、これを再利用することを考える。R からPython を利用するには、rPython パッケージを使う。
入らない場合は

sudo apt-get install python-dev

# これはPython
import numpy as np
from scipy.special import sph_harm as _sph_harm
def sph_harm(l, m, theta, phi):
    return _sph_harm(m, l, phi, theta)

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# R^2(u=θ, v=Φ)
u = np.linspace(0, np.pi, 80)
v = np.linspace(0, 2 * np.pi, 80)

# R^3
x = np.outer(np.sin(u), np.cos(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.cos(u), np.ones(v.size))

fig = plt.figure()
Lmax = 3
for l in range(Lmax + 1):
    for m in range(-l, l+1):
        # 球面調和関数の値を色で表す
        sph = np.abs(sph_harm(l, m, np.array(u, ndmin=2).T, np.array(v)))
# 正規化
#        if sph.max() == sph.min():
#            sph.fill(0.5)
#        else:
#            sph = (sph - sph.min()) / (sph.max() - sph.min())
        colors = cm.jet(sph) #cm.bwr
        ax = fig.add_subplot(Lmax + 1, 2*Lmax+1, l*(2*Lmax+1)+m+Lmax+1, projection="3d")
        ax.set_axis_off()
        surf = ax.plot_surface(x, y, z, facecolors=colors,
                rstride=2, cstride=2, linewidth=10, shade=False)

plt.savefig("sh_c.png", transparent=False)

rPython には4つの関数しかない。しかし、Python ができるのであればこれらの関数だけで十分強力な武器になる。
python.assign はPython へオブジェクトを定義する。

a <- 1:4
python.assign( "a", a )

python.exec はPython 上で演算を実行する。代入もできる。

python.exec( "b = len(a)" )

python.get はPython 上で定義されているオブジェクトを取り出す。ただし、Python でのarray 配列などはR に持ち込んだ時にJSON がウンタラカンタラ言われるので、Python 上でlist にしてR に吐き出すのがいいのかもしれない。

python.get("b")

python.call はPython 上の関数に、R から引数をいれて計算させる。python.assign でいちいちPython になげてやってもいいけれども、面倒だったらpython.exec でPython 上に好きに関数を定義してしまえばR から解決できる。

a <- 1:4
b <- 5:8
python.exec( "def concat(a,b): return a+b" )
python.call( "concat", a, b)

Python でライブラリを持っていれば、実行できる。

python.exec( "import math" )
python.get( "math.pi" )

さて、RでできないことをPython でゴリ押ししよう。scipy 内に球面調和関数を計算してくれるやつがあるので、これをそのまま読み込む。
リンクではlとm, $\theta$ と $\phi$ が逆なのでゴニョゴニョしているが、R で書く行が増えるのでいい感じに帳尻を合わせる。
Python で1行のものはそのままpython.* で挟めばよい。

# R で
# rPython
python.exec("import numpy as np")
python.exec("from scipy.special import sph_harm as _sph_harm")

Python はインデント構造にうるさいので、タブインデントをひたすら横に書いた長い行をいれてもよいが、見栄えを大切にするなら改行コマンドをいれる。
元のコードが

# Python
def sph_harm(l, m, theta, phi):
  sph = np.abs(_sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T))
  return [list(d0) for d0 in sph]

ゴリ押しすると

# R
python.exec("def sph_harm(l, m, theta, phi):\
  sph = np.abs(_sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T))\
  return [list(d0) for d0 in sph]")

また、上の実部だけ取り出すやつも、元のコードが

# Python
def sph_harm(l, m, theta, phi):
  if m > 0:
      sph = np.sqrt(2) * (-1)**m * _sph_harm( m, l, np.array(phi), np.array(theta, ndmin=2).T).real
  if m == 0:
      sph = _sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T).real
  else:
      sph = np.sqrt(2) * (-1)**m * _sph_harm(-m, l, np.array(phi), np.array(theta, ndmin=2).T).imag
  return [list(d0) for d0 in sph]

ゴリ押しすると

# R
python.exec("def sph_harm(l, m, theta, phi):\
  if m > 0:\
      sph = np.sqrt(2) * (-1)**m * _sph_harm( m, l, np.array(phi), np.array(theta, ndmin=2).T).real\
  if m == 0:\
      sph = _sph_harm(m, l, np.array(phi), np.array(theta, ndmin=2).T).real\
  else:\
      sph = np.sqrt(2) * (-1)**m * _sph_harm(-m, l, np.array(phi), np.array(theta, ndmin=2).T).imag\
  return [list(d0) for d0 in sph]")

これで、R上からm, l, $\theta$ , $\phi$ を与えればPython で球面調和関数を計算して、あたかもR でやってくれたようにできる関数ができたので、確かめる。

# R^2(u=θ, v=Φ)
# パラメータ
u <- seq(0, 180, length=180)*pi/180 # u = np.linspace(0, np.pi, 80)
v <- seq(0, 360, length=360)*pi/180 # v = np.linspace(0, 2 * np.pi, 60)

# プロット用の xyz
x <- c(outer(sin(u), cos(v)))
y <- c(outer(sin(u), sin(v)))
z <- c(outer(cos(u), rep(1, length(v))))
xyz <- cbind(x, y, z)

l <- 6
m <- 6

# ここで rPython の利用
# sph_harm は実数部分版
sph <- python.call("sph_harm", l, m, u, v) # u length list and v length vector
sph <- unlist(sph)
cuti <- 20
cols <- greenred(cuti)[cut(sph, quantile(sph, seq(0, 1, length=cuti + 1)))]

plot3d(xyz, col=cols)
zz <- t(polar2rect(abs(sph), t(expand.grid(v, u)[,2:1])))
plot3d(zz, col=cols)

　
ubuntu14.04, R3.3.0, Python 2.7.1 なので他OS では未確認。