numpy,scipy,sympyを実務に使って練習する - ryamadaのコンピュータ・数学メモ

ライブラリ・モジュール、その読み込み

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

- 頂点座標を持つ２次元アレイと頂点IDの三つ組を渡して、複数の三角形の重心を計算する(頂点IDトリオ、頂点座標２次元アレイのハンドリング関数は後述)

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

- 頂点座標を持つ２次元アレイと頂点IDの三つ組を渡して、複数の三角形の面積を計算する

# 三角形の面積を頂点ID三つ組リストと頂点座標格納行列とから作る
# 頂点IDは、三角形の個数だけ、長さ３の整数型リストを束ねた２次元アレイ
# 頂点座標は、頂点の個数の長さのベクトルを３つ束ねたもの
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個の要素を格納できる３次元アレイに格納するのが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()