Kumpei Shiraishi's blog

物理、プログラミング、クライミングに関する雑記

Pythonで球面調和関数を描画

目標

Pythonを使って、下図のようなものを描く。

Spherical Harmonics

結論

好きな自然数{l}{m}{m=-l,-l+1,\dots,l-1,1})を以下のコード冒頭で指定して、実行する。

Draw spherical harmonics with your preferred quant …

そして以下のような図を得る。

20170530_sph_harm_3_1

能書き

私は、この図を一枚描くのにかなりの時間を取られたので、バカな自分向けに、極めて初歩的な所から備忘録的にメモをしていく。

モジュールのimport

  • numpy:自明
  • scipy.special.sph_harm:球面調和関数。後述。
  • matplotlib.pyplot:自明
  • mpl_toolkits.mplot3d.Axes3D:3次元の図を描くのに必要のようだ。

scipysph_harmについて

公式のreferenceを見ると、引数は以下の4つを取る。これがWikipedia等で用いられているnotation(球面調和関数球面座標系と異なることが、混乱の原因その1だった。

  1. m:普通の{m}と一緒。即ち、{m=-l,\dots,l}となる整数。
  2. n:Wikipedia等では{l}とされている自然数。
  3. theta:定義域は[0,2*pi]。つまり、Wikipedia等では{\phi}とされている角度成分。
  4. phi:定義域は[0,pi]。つまり、Wikipedia等では{\theta}とされている角度成分。

Tipとしては、用意したデータセットの範囲をよく確認して、どっちの引数に食わせるか考えると良い。後に参考文献で挙げるページでは、再定義しているが、私にはこれは余計混乱するだけだった。

データセットの構成

サンプルコードでは、極座標の角度のデータ数を100とした。このことは、絵的に100個の点を用意したことにはならない。10,000個の点を用意しているのである。極座標変換の式

{\displaystyle
x = r\sin\theta\cos\phi \\
y = r\sin\theta\sin\phi \\
z = r\cos\theta
}

を見れば分かるように、{x}{y}{\theta}{\phi}の積で得られる。よって、絵的には10,000点を用意しているのである。{z}はそれに合わせて10,000点用意する。

これをコードにすると、直積を取るnp.outerを用いることになる。これにより、100行1列のベクトル{\sin\theta}と1行100列のベクトル{\cos\phi}の積である100行100列の行列が生成される({x}の場合)。

np.meshgridで鮮やかに書けるが、とりあえず1ステップずつ、自分が何をやっているのか分かる方が良いと思ってnp.outerを使った。

これに対応する形で球面調和関数も作らないといけない。まずthetaphiの直積も作り、それらをsph_harmに放り込んで、用意した。もちろん最初に指定した量子数を使った。

球面調和関数が何を表しているか

目的の「例の図」で、球面調和関数は何に対応しているのか?それは、中心からの距離である。半径rの3次元球に、生成した球面調和関数を掛けると、「例の図」の1パーツが出来上がる。

注意したいのは、一般に球面調和関数は複素数値を取るということだ。実部と虚部の物理的な意味はよく把握していないが、「例の図」を描くには実部を選べばよい。また、中心からの距離なので、実部を選んだあとで、絶対値を取る(順番を間違えないように)。

3次元描画

3次元描画の方法は、色々ある。まずfigureのインスタンスを生成、サブプロットでprojection='3d'とすれば、3次元で描ける。

サンプルで生成した図では、wireframeにした。他の事例は参考文献参照。

それを保存すれば、めでたしめでたし。

参考文献

matplotlibのデフォルトの色

2017年5月28日現在、matplotlibの最新版は2系である。このバージョンでは、デフォルトのプロットカラーが1系と変更されている。 2系でのデフォルトの色は以下の通りである。

番号 カラーコード サンプル
C0 #1f77b4
C1 #ff7f0e
C2 #2ca02c
C3 #d62728
C4 #9467bd
C5 #8c564b
C6 #e377c2
C7 #7f7f7f
C8 #bcbd22
C9 #17becf

例えば、

plt.plot(x, y, color='#ff7f0e')

と書けば、オレンジ色でグラフを描画できる。

'C0'から'C9'で色が指定できるなんて知らなかったが、trivialな番号を既知とする可読性の低いコードになると思うから、僕は使わないだろう)

ソース:Changes to the default style

DDSKKの統計をグラフにするコードを公開した

以前、sedとPython(pandas、matplotlib)を使って、DDSKKの統計をグラフにする手順をまとめた記事を書いた。

kumpeishiraishi.hatenablog.com

この記事に書いたような煩雑な手順を追わなくとも、Pythonとmatplotlibだけでグラフを描けるようにした。

github.com

詳しくはREADMEを参照。チンケなプログラムだが、あまり顧みられていない分野を自分で見つけてきて、それを埋めるプログラムを書いた初めての経験なので嬉しく思う(まあ、注目されていないということは、大して意味がないということではあるのだが)。皆さん、是非遊んでみてください。

SKKの統計をグラフにする

概要

SKKは、個人辞書に登録されている候補数の統計を記録することができる。この統計情報のグラフ化を試みる中で、以下のことを勉強することができた:

  • sedの使い方
  • pandasでcsvファイルを取り扱う方法
  • datetimeで日時データを取り扱う方法
  • matplotlib.datesで時系列データを表示する方法

この記事は、その個人的な記録である。グラフ化の元ネタはこのページ:~/.skk-record « 余談ですが……

結論

グラフ

Candidates in my SKK

コード・手順

  1. % cp ~/.skk-record ~/temp/rec.csvで、作業フォルダに統計情報をコピーする
  2. % sed -i '' -e 's/ 登録.*語数: /,/g' rec.csvで、日時と語数の列からなるcsvファイルに整形する(BSD sed)
  3. 以下のPythonコードを実行する:
import datetime as dt
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

# データフレーム
df = pd.read_csv('rec.csv',header=None)
df.columns = ['datetime','words']

# x軸時系列データのリストを整える
dates = [None]*df.shape[0]

for i in range(df.shape[0]):
    dstr = df.ix[i,'datetime']
    ddate = dt.datetime.strptime(dstr, '%a %b %d %H:%M:%S %Y')
    dates[i] = dt.datetime(ddate.year,ddate.month,ddate.day,ddate.hour,ddate.minute,ddate.second)
    # 当たり前だが、ここで時、分、秒を格納しなければ、もっと粗いグラフになる

# 以下、オブジェクト指向なmatplotlibの使い方をしている
# まず基本のFigureクラスのインスタンス。キャンバスだと思っておけばよい。
fig = plt.figure()

# 次にAxesクラスのインスタンス。白い背景と軸くらいに思っておこう。
ax = fig.add_subplot(1,1,1)

# axにプロットする。labelは凡例の名前を指定するため。
ax.plot(dates,df.ix[:,'words'],label='SKK jisyo candidates')

# 以下でグラフの体裁を整える
# 凡例の位置はlegend関数で指定できる
ax.legend(loc='upper left')

# axのx軸について、locatorとformatterを指定する
ax.xaxis.set_major_locator(mdates.DayLocator())
ax.xaxis.set_minor_locator(mdates.HourLocator(interval=12))
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d'))

# autofmt_xdate()で良い感じに回転させる
fig.autofmt_xdate(bottom=0.125,ha='center')

# 描画範囲を指定
start = min(dates)
end = max(dates)
ax.set_xlim(dt.date(start.year,start.month,start.day),dt.date(end.year,end.month,end.day+1))

# axにy軸のラベルを付加する
ax.set_ylabel('Number of words')

plt.savefig('a.png',dpi=300)

勉強したこと

sed

sedは置換、挿入、削除などのテキスト処理ができるコマンド。使い方は

を参照。ここでは今回使ったものだけを書いておく。

-eオプションは、その後に処理内容が来ることを表す。

-iオプションは、処理結果をファイルに上書きすることを表す。GNU sedとBSD sedで挙動が異なるので、後述する注意を参照。

特殊文字は、

  • .:任意の一文字
  • *:直前の文字が任意の個数続く文字列(0個も含む)
  • <space><space>*:連続するSPACEを表す
  • <tab>:TAB
  • .*:任意の文字列
  • \^:行頭
  • \$:行末

コマンドの前に来る数字をアドレスと呼び、その行のみに対して処理を行う。例:'3d'ならば3行目のみを削除する。

置換処理はs(Substitution)コマンドで、's/処理前文字列/処理後文字列/g'のように用いる。g(Global)オプションを渡さない場合、各行の最初に登場した当該文字列にのみ処理を行う。区切り文字はsコマンドの後に来たものなら何でも良い。

さて、sedにはGNU系とBSD系があり、それぞれで挙動が違うことに注意。Mac OS XにはBSD sedが入っているが、Homebrewでgnu-sedを入れればGNU sedを使うこともできる。BSD sedを使う上で注意が必要なのは、

  • -iオプションを付けて上書きする場合
  • 改行\nを出力する場合

である。

-iオプションに関してBSD sedのマニュアルには

Edit files in-place, saving backups with the specified extension.

と書いてある。つまりBSD sedにおいては、-iオプションの後ろで指定した文字列をファイル名の末尾に付けて元ファイルのバックアップを取るのである。例:sed -i 'boke' -e 's/foo/bar/g' hoge.txtを実行したならば、hoge.txtに置換済みの内容(欲しい内容)が保存され、hoge.txtbokeに置換前の内容(バックアップ)が保存される。GNU sedと同じように、余計なファイルを生成せずに上書きしたいなら-i ''と空の文字列を指定すればよい。

改行出力についてはこのページを参考にした。結論として、例えば1行目に文字列を挿入して上書きするには、

sed -i '' 1s/\^/foobar\\$'\n'/ input.txt

を実行しよう(今回はpandasでヘッダーを付加したので、BSD sedでの改行出力は不要なのだが)。

pandas

pandasはPythonでデータ分析をする際の定番ライブラリ。インストールしたものの、碌に使っていなかったので、ちょうど良い機会だった。

以下のページを主に参照した:

ここでは、今回のコードを読むのに必要なこと等を書いておく。

csvデータを読み込むには、read_csv('ファイル名')とする。第1行をヘッダーとして認識するので、これが不要ならheader=Noneを与える。また、読み込み時にヘッダーを指定するにはnames='ABC'を与える。そうすれば、第1列をAなどと名付けられる。

read_csv()でデータフレームdfを生成した後から、ヘッダーを指定して列名を名付ける(上書き)にはdf.columns = ['columnA','columnB','columnC']とする。

shape(行数, 列数)を、describe()で平均、分散などの基本的な統計量を得ることができる。

各要素にアクセスするには、ix[]を用いる。例:df.ix[2,'columnB']で第3行、第2列の要素を得る。列は列名でも列番号でも指定できる。複数選択は:を使う。

ソートはsort()。引数に列名を渡して、昇順でソートする。降順はascending=False

行の追加には、まず追加するデータフレームadd = DataFrame([['foobar',1]])を用意する。列名は追加先のデータフレームdfと同じものを、indexは追加先の末尾に合うものを指定する(然もなくばデータフレームの列が拡張され、形が変わる)。そしてconcat([df,add],axis=0)が追加後のデータフレームである。axis=1とすれば列の追加になる。

datetime

datetimeはPythonで日時データを扱うモジュール。使い方は

がよくまとまっている。

文字列から日時にするには

tstr = 'Sat Apr 23 07:12:25 2016'
tdatetime = dt.datetime.strptime(tstr, '%a %b %d %H:%M:%S %Y')

とする。日時から文字列にするには

tdatetime = dt.datetime(year=2015, month=10, day=21, hour=16, minute=29)
tstr = tdatetime.strftime('%Y/%m/%d')

とする。

指定子を知るにはここを参照。現在時刻を取得する.now()などの使い方も載っている。

matplotlib.dates

詳しくは以下を参照:

matplotlibは、オブジェクト指向なプロットを提供する。そのようにプロットする際の流れはコードのコメント、またはここに挙げたページを見てほしい。ここでは、matplotlib.datesの簡単な使い方について書いておく。

手作業で日時を目盛tickに指定するのはかなり大変だ。しかしmatplotlib.datesを使えば、tickの位置を制御するlocatorと、tick毎につける名前の書式を指定するformatterを、簡単に指定できる。

locatorでは分、時、日、週、月、年などの単位でインターバルを指定できる。例:matplotlib.dates.HourLocator(interval=6)は6時間毎にtickを記すlocatorである。他にも、毎週月曜日や、3年毎のイースターのようなtickを簡単に指定できる(公式ドキュメント参照)。

fotmatterでは、matplotlib.dates.DateFormatter()を使うことがほとんどだろう。例:matplotlib.dates.DateFormatter('%b-%d')はApr-24という形式で出力するformatterである。

これらをxaxis.set_major_locator()、xaxis.set_major_formatter()などに渡せばよい。

ちなみに、x軸のformatに日時を指定したが、多くの場合「2016/04/24」や「08:31:54」は横に長すぎて重なってしまう。autofmt_xdate()は、それを良い感じ回転、位置移動してくれる。デフォルトはbottom=0.2, rotation=30, ha='right'で、それぞれ下の余白、回転角、tickの相対位置を表す。また、今回はset_xlim()で描画範囲を指定した。指定しないと原点の下にx軸のtickが表示されず、ちょっと気持ち悪かったからだ。具体的には、x軸データの最大最小を取得し、両者の年月日のみを抜き出して、それを描画範囲とした。これならば、各日付の0時から24時まで描画されるので、欲しいtickになる。

matplotlibで定数関数を描く

やりたいこと

matplotlibで定数関数を描画したい。

ダメな例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 1

x = np.arange(-1,1,0.1)
y = f(x)
plt.plot(x,y)

plt.show()

これだと

ValueError: x and y must have same first dimension

というエラーが出る。xyのデータ数が合っていない、ということだ。

上手くいく例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 1

x = np.arange(-1,1,0.1)
y = [f(x[0])] * np.size(x)
plt.plot(x,y)
plt.show()

定数関数だと分かっているので、xの要素のどれかを指定してf(x)に渡し、xの配列長と同じ長さのリストをyに代入すれば良い。

ちなみにyの初期化は、こちらの方法の方がリスト内包表記よりも5.7倍くらい速かった。

matplotlibでグラフの場合分け

やりたいこと

matplotlibで、場合分けしたグラフを描きたい。

この記事では参考ページに倣って、

{
f(x) = \left\{ \begin{array}{ll}
x^2 & (0 \leq x \leq 1)\\
2-x & (1 < x \leq 2)\\
x^2 - 3x +2 & (2 < x \leq 3)\\
\end{array} \right.
}

を描く。

ダメな例

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    if 0<=x and x<=1:
        return x**2
    elif 1<x and x<=2:
        return 2-x
    elif 2<x and x<=3:
        return x**2-3*x+2

x = np.arange(0,3,0.001)
y = f(x)
plt.plot(x,y)

plt.show()

これを実行すると

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

というエラーになる。

関数f(x)の定義部分でxの大きさを判定しているのに、x = np.arange(0,3,0.001)というnumpyの配列を渡しているのが原因らしい。

上手くいく例

import numpy as np
import matplotlib.pyplot as plt

def f1(x):
    return x**2
def f2(x):
    return 2-x
def f3(x):
    return x**2-3*x+2

x1 = np.arange(0,1,0.001)
y1 = f1(x1)
plt.plot(x1,y1,"b")

x2 = np.arange(1,2,0.001)
y2 = f2(x2)
plt.plot(x2,y2,"b")

x3 = np.arange(2,3,0.001)
y3 = f3(x3)
plt.plot(x3,y3,"b")

plt.show()

明示的に別々の関数を定義すると上手くいく。複数の関数をmatplotlibで描画するには、単純にplotを複数回呼べばよい。グラフの色が別々だと格好悪いので、"b"を渡して全部青色にした。

そうして得られたグラフは以下。

Case in matplotlib

参考

3次元自己回避ランダムウォークをPythonで実装する

自己回避ランダムウォーク

自己回避ランダムウォーク(self avoiding random walk, SAW)とは、過去に占有した座標を回避しながら進むランダムウォークの一種である。3次元版のSAWは高分子のシミュレーションでよく用いられる。3次元では、以下のように漸近する。

{
\langle R^{2} (N) \rangle \sim N^{2\nu}
}

このとき{\nu \simeq 0.592}であることが知られている。

Rosenbluth-Rosenbluth最適化

さて、こんな3次元SAWをPythonでシミュレーションしてみようと思う。3次元直交座標を考え、原点から粒子がランダムウォークして動いていくとする。モノマー50個の長さのポリマーなら50回動くという要領で、動いた回数がポリマーの長さに対応するのである。

素朴な方法、次の格子をランダムに選び、それが自分の動いてきた履歴に含まれていたら終了、をそのまま実装したとする。ある程度の長さのポリマーを得るには、この方法では効率が悪すぎる。例えば長さ40のポリマーにまで到達する粒子は、わずか0.04%ほどである。そこで何らかの最適化が必要になる。

今回は、Rosenbluth-Rosenbluth最適化(Rosenbluth, Rosenbluth. J. Chem. Phys. 23, 356 (1955); http://dx.doi.org/10.1063/1.1741967)という方法を採る。この方法では次の格子をランダムに選ぶのではなく、6つの候補({x}正方向・負方向、{y}正方向・負方向、{z}正方向・負方向)から既に履歴に含まれているものを除外してランダムに選び、その代わりに重みをつける。

具体例を考えよう。{N-1}回の移動が終わり、{N}回目の移動を行う段階を考える。重み付けをする因子は{W(N)}である。6つの候補のうち、

  1. 履歴に含まれるものが1個(直前の位置)のみの場合、{W(N) = W(N-1)}
  2. 履歴に含まれるものが{m}個({m = 2, 3, 4, 5})の場合、{W(N) = \frac{m}{5}W(N-1)}
  3. 履歴に含まれるものが6個の場合、{W(N) = 0}として終了

以上の規則にしたがって、{W(N)}を決定する。但し、{W(0) = 1}である。

最後に、平均2乗変位{\langle R^{2} (N) \rangle}は以下の表式で求められる。

{ \displaystyle
\langle R^{2}(N) \rangle = \frac{\sum W(N)R^{2}(N)}{\sum W(N)}
}

但し、和は粒子番号についてとる。

Pythonコード

近日公開