yukiのブログ

作ったものなど

私のemacs環境をコピーする方法

emacsの設定をあれこれいじった結果それなりに使えるものになったので、emacsがインストールされていないPCで同じ設定のemacsを利用する方法をメモしておきます。未来の自分用です。

環境

Ubuntu 14.04 LTS

できるようになること

環境をコピーして拡張機能をすべて使えるようにすれば、C, OCaml, Verilog HDLでシンタックスハイライト、補完、リアルタイムなエラーチェックができるようになります。

コピーする手順

順に説明します。

emacsの最新版をインストー

Ubuntu 14.04 LTSでは普通にapt-getしても最新版のEmacsが入らないので、手動でppaを追加する必要があります。

sudo apt-add-repository -y ppa:adrozdoff/emacs
sudo apt update
sudo apt install emacs25

~/.emacsと~/.emacs.elを削除

そのままです。もしかしたらどちらも作られて無いかも。そうだったら何もしないでいいです。

リポジトリをクローンして、~/.emacs.d/におく

リポジトリここ。もともと存在するファイルは全て削除してください。

Caskをインストールする

Caskはemacs用のパッケージマネージャです。emacsとは別にインストールする必要があります。

インストール方法は公式サイトを参照。現在のバージョンでは以下のコマンドです。

curl -fsSL https://raw.githubusercontent.com/cask/cask/master/go | python

~/.emacs.d で cask install を実行

~/.emacs.d/Cask に使用するパッケージが羅列されており、cask install を実行することでそれらのパッケージがすべてインストールされます。

通常のインストールではcaskの実行ファイルにPATHが通らないっぽいので、その場合は

~/.emacs.d$ ~/.cask/bin/cask install

としてください。

emacsを起動して、エラーを潰していく

cask install が終わったあとにemacsを起動すれば基本的な拡張機能が読み込まれますが、一部は正常に動きません。一部の拡張機能には追加のパッケージが必要なためです(clang-formatterとかironyとかGNU Globalとか)。

具体的にどの拡張機能にどのパッケージが必要なのかは忘れたので、エラーを見てその都度パッケージ名で検索してください。そんなにアクロバティックなことはやってないので、ググって公式サイトのやり方に従えば大丈夫なはずです。

拡張機能の読み込みに関するエラーは*init log*タブか*minibuffer*タブに表示されます。Ctrl-x + Ctrl-b でタブの一覧を表示してから見に行ってください。

emacs拡張機能に必要となるパッケージは基本的にはapt-getで入るはずですが、Ubuntuのバージョンが最新でないために面倒な手順が必要になることがあります。特に GNU Global はソースからビルドする必要があるはずです。

多分この部分が一番めんどくさいでしょう。

(必要なら)フォントをインストールする

フォントにはrictyを使っています。インストールするには公式サイトからダウンロードして自分でビルドする必要があります。

別のフォントを使う場合には、~/.emacs.d/inits/98_themes.elを編集してください。

Sympyで方程式を解いてみる

Sympyで方程式を解いてみる

Sympyがまだ面白そうだったのでもっと色々やりました。

とりあえず前に受けた微分方程式講義ノートを引っ張り出してそこに書いてある練習問題を解いていきました。

微分方程式はすぐ解ける形もあればあれこれ手を焼いてやらないと解けないものもあるようです。すぐに解けないときは小手先の変更(hintを変えるとか)だけではうまくいかないことが多いので、さっさと地道な方法で解いたほうがいいと思います。

この記事はjupyter notebookで書いたのをmarkdownに変換したものです。notebookはここにあります。

# 初期処理
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
from  sympy.solvers.ode import *
init_printing()

# デフォルトだと微分方程式を解いたときに誤差関数(error function)を使うのでそれを書き直す関数
def expand_erfs(expr):
    """ 誤差関数を開く """
    return expr.replace(erf, lambda x: 2/sqrt(pi) * Integral(exp(-t**2), (t, 0, x)))\
               .replace(erfi, lambda x: 2/sqrt(pi) * Integral(exp(t**2), (t, 0, x))).simplify()

# 積分定数を潰す関数constantsimp()は積分定数を引数として
# 渡さなければいけないので、C1,C2,...をシンボルとして定義しなければいけない。
# それだとめんどくさいので、積分定数でないシンボルを渡せばすむ関数があると便利。
def constsimp(expr, variable):
    """ variable以外の自由変数をconstantsimp()に渡す """
    return constantsimp(expr, [c for c in ans.free_symbols if c != variable])

x = symbols("x")
a, b, c, t = symbols("a b c t")
y, f, g, P, Q, R = symbols('y f g P Q R', cls=Function)

方程式

expr = Eq(x**2+3*x-4)
display(expr, solveset(expr, x)) # solve()という関数もあるがsolveset()が推奨されている

\(\displaystyle x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-4, 1\right\}\)

expr = x**3+3*x**2+3*x-4
display(Eq(expr), solveset(expr))

\(\displaystyle x^{3} + 3 x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-1 + \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} - \frac{\sqrt{3} i}{2} \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} + \frac{\sqrt{3} i}{2} \sqrt[3]{5}\right\}\)

a, b, c = symbols("a b c")
solve(a*x**2+b*x+c, x)

\(\displaystyle \left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]\)

solveset(Gt(x**2, 0), domain=S.Reals)

\(\displaystyle \left(-\infty, 0\right) \cup \left(0, \infty\right)\)

solveset(sin(x) - sqrt(3)/2, domain=S.Reals)

\(\displaystyle \left\{2 n \pi + \frac{2 \pi}{3}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{3}\; \middle|\; n \in \mathbb{Z}\right\}\)

微分方程式

変数分離型

eq = Eq(Derivative(y(x)), g(x)*g(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = g{\left (x \right )} g{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \int^{y{\left (x \right )}} \frac{1}{g{\left (y \right )}}\, dy = C_{1} + \int g{\left (x \right )}\, dx\)

例題: ロジスティック方程式

# シンボルを定義
t, r, K = symbols("t r K")
N = symbols("N", cls=Function)
# ロジスティック方程式を定義
eq = Eq(Derivative(N(t)), r*(1-(N(t)/K))*N(t))
display(eq)
# 一般解を求める
ans = dsolve(eq)
display(ans)

# 特殊化する
args = [(K, 100), (r, 0.75)]
init_val = 1.0
func = ans.rhs.subs(args)
# 初期値を方程式の解として求める
c = solve(Eq(func.subs([(t, 0)]), init_val), domain=S.Reals)
# プロットする
C1 = symbols("C1") # 積分定数用のシンボルを定義しなきゃいけないのがめんどくさい
plot(func.subs([(C1, c[0])]), (t, -10, 15))

\(\displaystyle \frac{d}{d t} N{\left (t \right )} = r \left(1 - \frac{1}{K} N{\left (t \right )}\right) N{\left (t \right )}\)

\(\displaystyle N{\left (t \right )} = \frac{K e^{C_{1} K + r t}}{e^{C_{1} K + r t} - 1}\)

png

<sympy.plotting.plot.Plot at 0x23613929b70>

LiouVille(リウビル)型

LiouVille = Eq(y(x).diff(x,2) + P(x)*y(x).diff(x) + Q(y(x))*y(x).diff(x)**2)
display(LiouVille, dsolve(LiouVille, y(x), hint="Liouville"))

\(\displaystyle P{\left (x \right )} \frac{d}{d x} y{\left (x \right )} + Q{\left (y{\left (x \right )} \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int e^{- \int P{\left (x \right )}\, dx}\, dx + \int^{y{\left (x \right )}} e^{\int Q{\left (y \right )}\, dy}\, dy = 0\)

例題

eq = Eq(y(x).diff(x,2) + 2*x*y(x).diff(x) + 2*y(x)*y(x).diff(x) ** 2)
ans = dsolve(eq, hint="Liouville")
display(eq, ans, constsimp(expand_erfs(ans), x))

\(\displaystyle 2 x \frac{d}{d x} y{\left (x \right )} + 2 y{\left (x \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \operatorname{erf}{\left (x \right )} + \frac{\sqrt{\pi}}{2} \operatorname{erfi}{\left (y{\left (x \right )} \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int_{0}^{x} e^{- t^{2}}\, dt + \int_{0}^{y{\left (x \right )}} e^{t^{2}}\, dt = 0\)

eq = Eq(y(x).diff(x), sin(x)*tan(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \sin{\left (x \right )} \tan{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \left [ y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}\right ]\)

同次型

eq = Eq(y(x).diff(x), f(y(x)/x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = f{\left (\frac{1}{x} y{\left (x \right )} \right )}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- \int^{\frac{x}{y{\left (x \right )}}} \frac{f{\left (\frac{1}{u_{2}} \right )}}{u_{2} f{\left (\frac{1}{u_{2}} \right )} - 1}\, du_{2}}\)

eq = Eq(y(x).diff(x) , exp(y(x)/x) + y(x)/x)
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = e^{\frac{1}{x} y{\left (x \right )}} + \frac{1}{x} y{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \log{\left (\left(\frac{1}{C_{1} - \log{\left (x \right )}}\right)^{x} \right )}\)

分数型

eq = Eq(y(x).diff(x), (x-2*y(x)+3)/(2*x+y(x)-4))
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x - 2 y{\left (x \right )} + 3}{2 x + y{\left (x \right )} - 4}\)

\(\displaystyle \left [ y{\left (x \right )} = - 2 x - \sqrt{C_{1} + 5 x^{2} - 10 x} + 4, \quad y{\left (x \right )} = - 2 x + \sqrt{C_{1} + 5 x^{2} - 10 x} + 4\right ]\)

eq = Eq(y(x).diff(x), x*y(x)/(x**2+y(x)))
display(eq, dsolve(eq, hint="lie_group"))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x y{\left (x \right )}}{x^{2} + y{\left (x \right )}}\)

\(\displaystyle \left [ y{\left (x \right )} = C_{1} \left(- \sqrt{- 2 C_{1} x^{2} + 1} - 1\right), \quad y{\left (x \right )} = C_{1} \left(\sqrt{- 2 C_{1} x^{2} + 1} - 1\right)\right ]\)

# 自動で解けないのもあるので、そのときは答を代入する
eq = Eq(y(x).diff(x), x**2*y(x)/(x**3+y(x)))
# この答えは手で解いたもの
ans_implicit = Eq(y(x)**3/(2*x**3+3*y(x)), 6*C1)
# 実数解だけ取る(虚数解でも成り立つだろうけどとりあえず無視)
ans = solve(ans_implicit, y(x))[0]
# これがTrueなのでans_implicitが正しいと分かる
display(eq, ans_implicit, Eq(y(x), ans), eq.subs([(y(x), ans)]).doit().simplify())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x^{2} y{\left (x \right )}}{x^{3} + y{\left (x \right )}}\)

\(\displaystyle \frac{y^{3}{\left (x \right )}}{2 x^{3} + 3 y{\left (x \right )}} = 6 C_{1}\)

\(\displaystyle y{\left (x \right )} = - \frac{18 C_{1}}{\sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}} - \frac{1}{3} \sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}\)

\(\displaystyle \mathrm{True}\)

定数変化法

eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq, dsolve(eq, y(x), hint="almost_linear"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \int Q{\left (x \right )} e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \int P{\left (x \right )}\, dx}\)

eq = Eq(x*y(x).diff(x) - 2*y(x), x**3*cos(x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle x \frac{d}{d x} y{\left (x \right )} - 2 y{\left (x \right )} = x^{3} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

# 一般解から解いてみる
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq.subs([(P(x), -2/x), (Q(x), x**2*cos(x))]),
        dsolve(eq, hint="almost_linear").subs([(P(x), -2/x), (Q(x), x**2*cos(x))]).doit())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} - \frac{2}{x} y{\left (x \right )} = x^{2} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

Bernoulli(ベルヌーイ)型

n = symbols("n")
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x)*y(x)**n)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )} y^{n}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(\left(C_{1} - \left(n - 1\right) \int Q{\left (x \right )} \left(e^{- n \int P{\left (x \right )}\, dx}\right) e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \left(- n + 1\right) \int P{\left (x \right )}\, dx}\right)^{\frac{1}{- n + 1}}\)

eq = Eq(y(x).diff(x)+x*y(x), x*y(x)**3)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle x y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = x y^{3}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \frac{1}{\sqrt{\left(C_{1} + e^{- x^{2}}\right) e^{x^{2}}}}\)

Riccati(リッカチ)型

# 自動で解くのは無理っぽい
eq = Eq(y(x).diff(x) - y(x)**2 - (2-x)*y(x) + 2*x - 1)
spec_ans = x
display(eq.subs([(y(x), spec_ans)]).doit().simplify())
bern = eq.subs([(y(x), P(x)+spec_ans)]).doit().simplify()
ans = dsolve(bern, hint="Bernoulli")
display(eq, bern, ans, ans.subs([(P(x), y(x)-x)]))

\(\displaystyle \mathrm{True}\)

\(\displaystyle 2 x - \left(- x + 2\right) y{\left (x \right )} - y^{2}{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} - 1 = 0\)

\(\displaystyle - x P{\left (x \right )} - P^{2}{\left (x \right )} - 2 P{\left (x \right )} + \frac{d}{d x} P{\left (x \right )} = 0\)

\(\displaystyle P{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

\(\displaystyle - x + y{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

Clairaur(クレロー)型

eq = Eq(y(x), x*y(x).diff(x) - y(x).diff(x)**2)
dsolve(eq)

\(\displaystyle y{\left (x \right )} = \frac{x^{2}}{4} - \frac{1}{4} \left(C_{1} + x\right)^{2}\)

減衰振動

m, k, u = symbols("m k u")
eq = Eq(m*y(x).diff(x,2), -k*y(x) - u * y(x).diff(x))
dsolve(eq)

\(\displaystyle y{\left (x \right )} = C_{1} e^{\frac{x}{2 m} \left(- u - \sqrt{- 4 k m + u^{2}}\right)} + C_{2} e^{\frac{x}{2 m} \left(- u + \sqrt{- 4 k m + u^{2}}\right)}\)

定数係数斉次線形

# 一番簡単なやつ
eq = Eq(y(x).diff(x,2) - y(x).diff(x) - 6*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 6 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{3 x}\)

# 重解のパターン
eq = Eq(y(x).diff(x, 2) - 6*y(x).diff(x) + 9*y(x))
display(eq, dsolve(eq))

\(\displaystyle 9 y{\left (x \right )} - 6 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{3 x}\)

eq = Eq(y(x).diff(x,3) - 5*y(x).diff(x,2) + 8*y(x).diff(x) - 4*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 4 y{\left (x \right )} + 8 \frac{d}{d x} y{\left (x \right )} - 5 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{3}}{d x^{3}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \left(C_{2} + C_{3} x\right) e^{x}\right) e^{x}\)

eq = Eq(y(x).diff(x,4) - 8*y(x).diff(x,2) + 16*y(x))
display(eq, dsolve(eq))

\(\displaystyle 16 y{\left (x \right )} - 8 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{4}}{d x^{4}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- 2 x} + \left(C_{3} + C_{4} x\right) e^{2 x}\)

定数係数非斉次線形

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), exp(3*x))
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{3 x}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{e^{3 x}}{20}\)

eq = Eq(y(x).diff(x, 2) + 2*y(x).diff(x) + y(x), exp(-x))
display(eq, dsolve(eq))

\(\displaystyle y{\left (x \right )} + 2 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{- x}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x + \frac{x^{2}}{2}\right) e^{- x}\)

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) - y(x).diff(x) - 2*y(x), x**2)
display(eq, dsolve(eq))

\(\displaystyle - 2 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x^{2}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- x} + C_{2} e^{2 x} - \frac{x^{2}}{2} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) + 2*a*y(x).diff(x) + a**2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle a^{2} y{\left (x \right )} + 2 a \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- a x} + \frac{x}{a^{2}} - \frac{2}{a^{3}}\)

Sympyを試す

 

Sympyを試す

Sympyが面白かったので色々試しました。

このページはjupyter notebookを作ったものをmarkdownに変換したものです。もととなったnotebookはここにあります。

初期設定

# ここはとりあえずコピペすればよさそう
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
init_printing()

シンボルの定義

x = Symbol("x")
y = Symbol("y")            # これでもいいけど...
x, y, z = symbols("x y z") # こっちのほうが楽なのでこっちでいいと思う
f, g = symbols('f g', cls=Function)

# displayはnotebook側の関数で、引数を良い感じに出力してくれる。
display(x, y, z)

\(\displaystyle x\)

\(\displaystyle y\)

\(\displaystyle z\)

式の定義

expr = (x + y) ** 5

# display()を使うとnotebookを通じて出力できる
display(expr)
display(x * expr * y) # 変数の並び順は簡略化されるらしい

\(\displaystyle \left(x + y\right)^{5}\)

\(\displaystyle x y \left(x + y\right)^{5}\)

式の展開/因数分解

expr = (x + y) ** 5

display(expr.expand())                                # expand()で式を展開できる
display((x**3 + 3*x**2*y + 3*x*y**2 + y**3).factor()) # factor()はexpand()の逆

\(\displaystyle x^{5} + 5 x^{4} y + 10 x^{3} y^{2} + 10 x^{2} y^{3} + 5 x y^{4} + y^{5}\)

\(\displaystyle \left(x + y\right)^{3}\)

式の比較

expr = (x + y) ** 2
expanded = expand(expr)
display(expr)
display(expanded)

\(\displaystyle \left(x + y\right)^{2}\)

\(\displaystyle x^{2} + 2 x y + y^{2}\)

display(expr == expanded)              # この比較は代数的な比較ではない
display(Eq(expr, expanded))            # Eq()を使えば代数的な比較を作れる
display(Eq(expr, expanded).simplify()) # さらにSimplify()で恒等的に成り立つかどうかをチェックできる
False

\(\displaystyle \left(x + y\right)^{2} = x^{2} + 2 x y + y^{2}\)

\(\displaystyle \mathrm{True}\)

expr.equals(expanded) # equals()はランダムな点でテストをして等しいかどうかを判定するらしい
True

代入

expr = sin(x) + 1
expr.subs(x, pi / 3)

\(\displaystyle \frac{\sqrt{3}}{2} + 1\)

expr = x**2 + 3*x + 4
expr.subs(x, expr) # 式の代入も可能

\(\displaystyle 3 x^{2} + 9 x + \left(x^{2} + 3 x + 4\right)^{2} + 16\)

expr = x**3 + 5*x*y + x*z + 5*y**2 + y*z**2 + z**2
display(expr)
display(expr.subs([(x, 1), [y, 3], [z, 2]]))
display(expr.subs([(x, x), [y, x], [z, x]]))

\(\displaystyle x^{3} + 5 x y + x z + 5 y^{2} + y z^{2} + z^{2}\)

\(\displaystyle 79\)

\(\displaystyle 2 x^{3} + 12 x^{2}\)

無理数の計算

pi.evalf(30) # 1000ぐらい余裕

\(\displaystyle 3.14159265358979323846264338328\)

expr = sin(1)**2 + cos(1)**2
display((expr - 1).evalf())
display((expr - 1).evalf(chop = True))

\(\displaystyle -4.0 \cdot 10^{-124}\)

\(\displaystyle 0\)

numpy関数の利用

import numpy 
a = numpy.arange(100) 
expr = sin(x)**2 + cos(x)**3
display(expr)
f = lambdify(x, expr, "numpy") # numpyを使わないと遅い
f(a)

\(\displaystyle \sin^{2}{\left (x \right )} + \cos^{3}{\left (x \right )}\)

array([ 1.        ,  0.86580202,  0.75475425, -0.95036208,  0.29348079,
        0.94236043,  0.96327991,  0.86012577,  0.97574947, -0.58654075,
       -0.29478182,  0.9999805 ,  0.88881004,  0.92378616,  0.98385952,
       -0.01556116, -0.79539238,  0.90345118,  0.85189192,  0.98895835,
        0.90142733,  0.53566986, -0.99980414,  0.5648118 ,  0.89639378,
        0.99135691,  0.85223411,  0.88972231, -0.81857033,  0.02180437,
        0.97987667,  0.92866036,  0.88463129,  0.99982139, -0.33110279,
       -0.55466852,  0.98152993,  0.86256579,  0.95901972,  0.94785945,
        0.25853532, -0.93733553,  0.7760188 ,  0.86290782,  0.99984336,
        0.86900634,  0.73250099, -0.96191188,  0.32789382,  0.93680443,
        0.96737779,  0.85798037,  0.96910399, -0.61757906, -0.25811662,
        0.99952124,  0.89314651,  0.91891613,  0.98748892, -0.0530218 ,
       -0.77101524,  0.91618973,  0.85189896,  0.98629159,  0.90661841,
        0.50571095, -0.99823778,  0.59310782,  0.89154121,  0.99347747,
        0.85292682,  0.87499481, -0.84050791,  0.05901876,  0.97557656,
        0.93351789,  0.88062925,  0.99901083, -0.36702009, -0.5220171 ,
        0.98646955,  0.8652878 ,  0.95461493,  0.95326866,  0.22310345,
       -0.92285576,  0.79628569,  0.86033431,  0.99937407,  0.87250858,
        0.70927031, -0.97196403,  0.3617301 ,  0.93122345,  0.97129634,
        0.85614111,  0.96157091, -0.64772986, -0.22116672,  0.99847744])

極限/微分/積分

expr = sin(x)/x
# 極限をその場で計算しないならLimitを、するならlimitを使う(以下同様に2種類の関数がある)
display(Eq(Limit(expr, x, 0), limit(expr, x, 0)))
# 極限の方向は片方しか指定できないみたい
display(Eq(Limit(expr, x, 0, dir="-"), limit(expr, x, 0, dir="-"))) 

\(\displaystyle \lim_{x \to 0^+}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)

\(\displaystyle \lim_{x \to 0^-}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)

expr = exp(sin(x))
Eq(Derivative(expr), diff(expr)) 

\(\displaystyle \frac{d}{d x} e^{\sin{\left (x \right )}} = e^{\sin{\left (x \right )}} \cos{\left (x \right )}\)

expr = x * (exp(x)*sin(x) + exp(x)*cos(x))
Eq(Integral(expr), integrate(expr))

\(\displaystyle \int x \left(e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}\right)\, dx = x e^{x} \sin{\left (x \right )} - \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}\)

expr = exp(-x**2)
# 無限大は"oo"で指定
Eq(Integral(expr, (x, -oo, oo)), integrate(expr, (x, -oo, oo)))

\(\displaystyle \int_{-\infty}^{\infty} e^{- x^{2}}\, dx = \sqrt{\pi}\)

級数

expr = exp(x)
Eq(expr, expr.series(x=x, x0=0, n=10))

\(\displaystyle e^{x} = 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \frac{x^{6}}{720} + \frac{x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{x^{9}}{362880} + \mathcal{O}\left(x^{10}\right)\)

expr = exp(sin(x)**2)
Eq(expr, expr.series(x=x, x0=0, n=10))

\(\displaystyle e^{\sin^{2}{\left (x \right )}} = 1 + x^{2} + \frac{x^{4}}{6} - \frac{11 x^{6}}{90} - \frac{71 x^{8}}{2520} + \mathcal{O}\left(x^{10}\right)\)

プロット

x = symbols('x')
plot(x*x)
plot(sin(x))

png

png

<sympy.plotting.plot.Plot at 0x271b9e4d828>

正規表現の実装(python)

前回の記事で作ったεNFAを使って正規表現を実装できたのでまとめます。ソースは前回と同じくこのリポジトリに。

実装した正規表現

正規表現のすべての機能を実装するのは大変なので、実装するのは基本的なものだけにしました。

実装したのは

  • 一文字 (“a"など)
  • 文字の連結 (“ab"など)
  • 文字の和集合 (“(a|b)"など)
  • 表現の繰り返し(“a*"など)
  • 任意の一文字 (“."のこと)

のみです。教科書には定義に"Φ"(空文字列)がありましたが使い所がなさそうなので代わりに".“を実装しました。

“a+”(aの一回以上の繰り返し)や"(^a|b)“(aでもbでもない文字全て)などの機能は上の表現を組み合わせれば作れます。

他の特徴として入力される文字集合を指定しなければならないというのがあります。これは正規表現の都合というよりもこれまで作ってきたεNFAが文字集合がないと定義できないという都合のためです。

インターフェース

関数regex_to_eNFA()で正規表現をεNFAに変換し、eNFA.run()で文字列が正規表現にマッチするかどうかを調べます。

enfa = regex_to_eNFA("10*1", {"0", "1"})
enfa.run("101")   # -> True
enfa.run("10101") # -> False
enfa.run("10001") # -> True

実装の手順

正規表現をεNFAに変換するプログラムは大きく分けて3つの部分からなります。

  • 正規表現を一つづつに分ける
  • 一つづつに分けた正規表現をεNFAに変換する
  • 変換されたεNFAをε遷移で繋いで一つのεNFAを作る

この部分をコードにすると

def regex_to_eNFA(string, alphabet):
    """ 正規表現stringを認識するeNFAを返す """
    return eNFA.serial_connect([single_regex_to_eNFA(single_regex, alphabet) for single_regex in concat_split(string)])

となります。一行にまとめているので読みにくいですが、concat_split()で正規表現を一つづに分け、single_regex_to_eNFA()で一つづつに分けた正規表現をεNFAに変換し、最後にeNFA.serial_connect()でεNFAを繋いでいます。

順に説明します。

concat_split()

concat_split()は正規表現を一つづつに分ける関数です。一つづつと言うのは連結されていないという意味です。例えば、"a(b|c)de“は["a”, “(a|b)”, “d”, “e”]と分かれます。

ここでの詰まったのは、処理を始める前に一番外側のはずせるカッコを外しておく必要がある点です。もしカッコを外さないと、例えば"(a(b|c)d)“という正規表現は[”(a(b|c)d)“]としか分解されず、無限ループに陥ってしまうからです。対策としてremove_outer_branket()を作り、これに通してから処理を始めるようにしました。

def remove_outer_branket(string):
    """ 一番外側のカッコが外せるなら外して、外せないならそのまま返す """
    if not (string[0] == "(" and string[-1] == ")"):
        return string

    blanket_count = 0
    for char in string[:-1]:
        if char == "(":
            blanket_count += 1
        elif char == ")":
            blanket_count -= 1
        if blanket_count == 0:
            return string
        if blanket_count == 1 and char == "|":
            return string
    return remove_outer_branket(string[1:-1])

def concat_split(string):
    """ 正規表現の連結stringを正規表現一つづつに切って返す """
    ans = []
    right = left = 0
    string = remove_outer_branket(string)
    while left < len(string):
        if string[left] == "(":
            blanket_count = 1
            while blanket_count > 0:
                left += 1
                if string[left] == "(":
                    blanket_count += 1
                elif string[left] == ")":
                    blanket_count -= 1
        if left + 1 < len(string) and string[left + 1] == "*":
            left += 1
        ans.append(string[right:left + 1])
        right, left = left + 1, left + 1
    return ans

テストしてみます。なおgituhubのソースではdoctestを使ってテストをコメントに書いています。改変したときにあってるかどうかがすぐに調べられて便利です。

>>> remove_outer_branket("(abc)")
'abc'
>>> remove_outer_branket("(((abc))")
'abc'
>>> remove_outer_branket("(a|b)(c|d)")
'(a|b)(c|d)'
>>> concat_split("ab")
['a', 'b']
>>> concat_split("a*")
['a*']
>>> concat_split("(abc)*(a|b|c)(cba)*abc")
['(abc)*', '(a|b|c)', '(cba)*', 'a', 'b', 'c']
>>> concat_split("...")
['.', '.', '.']
>>> concat_split("((a|A)(b|B)(c|C))")
['(a|A)', '(b|B)', '(c|C)']
>>> concat_split("((a|A)(b|B)(c|C))*")
['((a|A)(b|B)(c|C))*']

single_regex_to_eNFA()

single_regex_to_eNFA()は連結のない正規表現を受け取り、それを認識するεNFAを返す関数です。

連結のない正規表現

  • 一文字 (“a"など)
  • 文字の和集合 (“(a|b)"など)
  • 表現の繰り返し(“a*"など)
  • 任意の一文字 (“."のこと)

のどれかなので、場合分けして処理します。

def is_dot(string):
    """ stringが任意の一文字を表す記号"."ならTrue """
    return string == "."

def is_atom(string):
    """ stringが通常の文字一つならTrue、それ以外はFalse """
    return True if len(string) == 1 and string not in "()|*." else False

def is_repeat(string):
    """ stringが"X*"ならTrue、それ以外はFalse """
    return len(string) > 1 and string[-1] == "*"

def is_union(string):
    """ stringが文字の集合を表す正規表現ならTrue """
    return string[0] == "(" and string[-1] == ")"

def single_regex_to_eNFA(string, alphabet):
    """ 一項からなる正規表現stringを認識するeNFAを返す """
    if is_dot(string):
        return eNFA.any_word(alphabet)
    if is_atom(string):
        return eNFA.one_word(string, alphabet)
    elif is_repeat(string):
        return eNFA.repeat(regex_to_eNFA(string[:-1], alphabet))
    elif is_union(string):
        return eNFA.parallel_connect([regex_to_eNFA(sub_regex, alphabet) for sub_regex in union_split(string)])
    else:
        assert False, "invalid regular expression: %s" % string

わかりやすいようにis_…()という関数を作っていますが、判定は文脈に依存しています。例えばis_union()は場合分けの最後に置かれなければ正しく動きません。なので関数にしなくても良かったかもしれません。

εNFAを作る関数

any_word(), one_word(), repeat(), parallel_connect(), serial_connact() は全てεNFAを作る関数です。一例としてserial_connect()を示し、解説は省略します。

def serial_connect(automaton):
    """ 複数のオートマトンを横並びに一つにまとめる """
    alphabet = set()
    states = set()
    transitions = {}
    for i, automata in enumerate(automaton):
        # alphabetを更新
        alphabet = alphabet.union(automata.alphabet)

        # statesを更新
        # 状態とautomataのインデックスを組にすることで唯一性を確保
        states = states.union(set([(i, state) for state in automata.states]))

        # automata内部の遷移をtransitionsに追加
        for final_state in automata.final_states:
            transitions[(i, final_state)] = {}
        for state, trans_dict in automata.transitions.items():
            transitions[(i, state)] = {}
            for char, states_transit_to in trans_dict.items():
                transitions[(i, state)][char] = frozenset([(i, state) for state in states_transit_to])

        # 一つ前のautomataからのイプシロン遷移を追加
        if i == 0:
            continue
        for pre_final_state in automaton[i - 1].final_states:
            transitions[(i - 1, pre_final_state)][-1] = frozenset([(i, automata.init_state)])

    # 始状態と終状態を作る
    init_state = (0, automaton[0].init_state)
    final_states = set([(len(automaton) - 1, f) for f in automaton[-1].final_states])
    return NFAWithEpsilonTransition(states, alphabet, transitions, init_state, final_states)

DFAへの変換と最小化

εNFAのDFAへの変換とDFAの最小化はこれまでに作ったので、正規表現から作ったεNFAを最小DFAに変換することができます。

これによって正規表現の性質が分かります。

>>> from Regex import regex_to_eNFA
>>> any = regex_to_eNFA("(0|1)*", {"0", "1"})
>>> any_mindfa = any.convert_to_DFA().minimized()
Automata
    states      : {{{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}
    alphabet   : {'0', '1'}
    transitions : {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}} : {'0': {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}, '1': {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}
    init_state  : {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}
    final_states: {{{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}

とても見にくいですが、よく見ればこのDFAは状態が一つしか無い、どんな文字列も認識するDFAだとわかります。正規表現"(0|1)“が表すのは任意の文字列なので、このDFAが正しく変換されていることがわかります。

最小DFAは全て同型なので、DFAが同型かどうかを判定する関数を作れば正規表現が同値かどうかを判定できるようになるはずです(やってないけど)。

またDFAの表記がみにくい問題については、DFAのコンストラクタで状態の名前を振り直す処理を加えればマシになるはずです。デバッグがしにくくなるので実装はしませんでしたが。

オートマトンをpythonで書く(続き)

この前作ったDFAをいくらか拡張したのでまとめます。ソースはgithubに。

NFA

講義ではDFAの次に非決定性有限オートマトン(Nondeterministic Finite Automata, NFA)が出てきたので、最初の拡張としてNFAを作りました。

NFAは文字列を読み込んでいる途中に状態を複数取ることができます。状態を複数取っているときに文字を読んだ場合、次の状態は現在の各状態の遷移先の和集合です。また文字列を認識するかどうかは最後まで文字を読み込んだあとの状態に終状態が一つでもあるかどうかで判断します。

こう書くと難しく感じますが、DFAのプログラムと骨格はほとんど同じなので特に迷わずに書くことができました。

class NondeterministicFiniteAutomata(Automata):
    """ 非決定性有限オートマトン """

    def run(self, string):
        """ stringを認識するかチェックする """
        states = {self.init_state}
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                states = self.next_states(states, char)
        return len(states.intersection(self.final_states)) != 0

    def next_states(self, states, char):
        """ statesからcharを読んだ時の次の状態集合を返す """
        next_states = set()
        for state in states:
            next_states = next_states.union(self.transitions[state].get(char, frozenset()))
        return frozenset(next_states)

参考までに、前の記事で書いたDFAのrun関数も載せておきます。NFAのrun関数と同じような形をしているのが分かるはずです。

class DeterministicFiniteAutomata(Automata):
    """ 決定性有限オートマトン """

    def run(self, string):
        """ stringを認識するかチェックする """
        state = self.init_state
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                state = self.transitions[state][char]
        return state in self.final_states

具体的なDFAは次のように作ります。

from Automata import NondeterministicFiniteAutomata as NFA
# 010を部分列に持つ文字列を認識するNFA
a, b, c, d = range(4)
states = {a, b, c, d}
alphabets = {"0", "1"}
transitions = {
    a: {"0": {a, b}, "1": {a}},
    b: {"1": {c}},
    c: {"0": {d}},
    d: {"0": {d}, "1": {d}},
}
init_state = a
final_states = {d}
has_010 = NFA(states, alphabets, transitions, init_state, final_states)

作った上でrun()すれば動作を確認できます。

# pythonコンソール
>>> has_010.run("0001011")
True
>>> has_010.run("0011")
False
>>> has_010.run("001100")
False
>>> has_010.run("100")
False
>>> has_010.run("0100")
True

NFAからDFAへの変換

任意のNFAを等価なDFAに変換するアルゴリズムがあるのでそれも作りました。状態の冪集合を考えることで、NFAが文字を読む間に取りうる状態の集合をDFAの一つの状態としてしまうことがこのアルゴリズムのポイントです。

下のコードでは目標となるDFAのstates, alphabets, transitions, init_state, final_statesを作ることでDFAを作っています。これらの中でalphabetsとinit_stateは自明であり、final_statesはstatesがわかってしまえばすぐに分かるので、作るときに苦労するのはstatesとtransitionsです。

    def convert_to_DFA(self):
        """ 等価なDFAに変換する """
        # 状態と遷移関数を作る
        states = {frozenset({self.init_state})}
        transitions = {}
        to_search = {frozenset({self.init_state})}
        searched = set()
        while len(to_search) != 0:
            searching = to_search.pop()
            searched.add(searching)
            for char in self.alphabets:
                # 実際に遷移させて考える
                reachables = self.next_states(searching, char)
                # 遷移先の状態集合そのものが一つの状態となる
                states.add(reachables)
                if reachables not in searched:
                    # 遷移させた先がまだ調べていない状態だったらさらに調べる必要がある
                    to_search.add(reachables)
                if not transitions.get(searching):
                    # 二重辞書の要素を一気に作ることはできないので注意
                    transitions[searching] = {}
                # 遷移を追加する
                transitions[searching][char] = frozenset(reachables)

        # DFAの終状態は終状態を一つでも含む状態全体からなる集合
        final_states = set([x for x in states if len(self.final_states.intersection(x)) != 0])
        # DFAの始状態はNFAの始状態だけからなる集合
        init_state = frozenset({self.init_state})
        # DFAのアルファベットはNFAと同じ
        alphabets = self.alphabets
        return DeterministicFiniteAutomata(states, alphabets, transitions, init_state, final_states)

次のように使います。同じ機能を持つDFAができていますが、変換によって状態数が増え、複雑になっていることがわかります。

# has_010を作った状態のpythonコンソール
>>> print(has_010)

Automata
    states      : {0, 1, 2, 3}
    alphabets   : {'1', '0'}
    transitions : 0 : {'1': {0}, '0': {0, 1}}
                  1 : {'1': {2}}
                  2 : {'0': {3}}
                  3 : {'1': {3}, '0': {3}}
    init_state  : 0
    final_states: {3}
>>> has_010_DFA = has_010.convert_to_DFA()
>>> print(has_010_DFA)
Automata
    states      : {{0, 2, 3}, {0, 2}, {0, 1, 3}, {0}, {0, 3}, {0, 1}}
    alphabets   : {'1', '0'}
    transitions : {0, 2, 3} : {'1': {0, 3}, '0': {0, 1, 3}}
                  {0, 2} : {'1': {0}, '0': {0, 1, 3}}
                  {0} : {'1': {0}, '0': {0, 1}}
                  {0, 1, 3} : {'1': {0, 2, 3}, '0': {0, 1, 3}}
                  {0, 3} : {'1': {0, 3}, '0': {0, 1, 3}}
                  {0, 1} : {'1': {0, 2}, '0': {0, 1}}
    init_state  : {0}
    final_states: {{0, 2, 3}, {0, 3}, {0, 1, 3}}

>>> has_010_DFA.run("010")
True
>>> has_010_DFA.run("0110")
False
>>> has_010_DFA.run("100110")
False
>>> has_010_DFA.run("1001010")
True

全てのNFAはDFAに変換できることと全てのDFAはNFAなことから、DFAとNFAが同じ計算能力を持つことがわかります。

εNFA

NFAをさらに拡張してイプシロン動作付き非決定性有限オートマトン(εNFA)も作りました。

εNFAとNFAの違いは、イプシロン遷移があるかどうかです。イプシロン遷移が状態\(p\)から状態\(q\)に張られていた場合、状態\(p\)に到達したら無条件で状態\(q\)にも到達したとみなします。さらに状態\(q\)から別の状態にイプシロン遷移が張られていれば、その状態についても同様です。

違いがこれだけなのでNFAをεNFAにするのは簡単かと思いましたがそれなりに手こずりました。イプシロン遷移を一回たどるだけでは不十分で、飽和するまでたどりきらないといけないためです。

class NFAWithEpsilonTransition(Automata):
    """ ε動作付き非決定性有限オートマトン """

    def reachables_with_a_epsilon_from(self, state):
        """ stateから一回のε遷移だけで到達可能な状態の集合を返す """
        return frozenset(self.transitions[state].get(-1, {}))

    def reachables_with_epsilons_from(self, state):
        """ stateからε遷移だけで到達可能な状態の集合を返す """
        ans = frozenset([state])
        reachables = self.reachables_with_a_epsilon_from(state)
        while not reachables.issubset(ans): # 飽和するまで繰り返す
            ans = ans.union(reachables)
            new_reachables = frozenset()
            for reachable in reachables:
                new_reachables = new_reachables.union(self.reachables_with_a_epsilon_from(reachable))
            reachables = new_reachables
        return ans

    def next_states(self, states, char):
        """ statesからcharを読んだ時の次の状態集合を返す """
        next_states = frozenset()
        for state in states:
            for reachable in self.transitions[state].get(char, frozenset()):
                next_states = next_states.union(self.reachables_with_epsilons_from(reachable))
        return frozenset(next_states)

    def run(self, string):
        """ stringを認識するかチェックする """
        states = self.reachables_with_epsilons_from(self.init_state)
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                states = self.next_states(states, char)
        return len(states.intersection(self.final_states)) != 0

以下のようにして動かします。

>>> from Automata import NFAWithEpsilonTransition as eNFA
>>> a, b, c, d, e = range(5)
>>> states = {a, b, c, d, e}
>>> alphabets = {"0", "1", "2", "3", "4"}
>>> transitions = {
...     a: {"0": {a}, -1: {b}},
...     b: {"1": {b}, -1: {c}},
...     c: {"2": {c}, -1: {d}},
...     d: {"3": {d}, -1: {e}},
...     e: {"4": {e}}
>>> }
>>> init_state = a
>>> final_states = {e}
>>> is_increasing = eNFA(states, alphabets, transitions, init_state, final_states)

>>> is_increasing.run("01144")
True
>>> is_increasing.run("01231")
False
>>> is_increasing.run("00004")
True
>>> is_increasing.run("33444")
True

NFAと同様の方法でεNFAもDFAに変換することができます(コードはほとんど同じなので省略します)。

DFAの最小化

任意のDFAは最小化することができます。つまり、任意のDFAに対して、そのDFAが同じ言語を認識するDFAのうち状態数が最小のDFAを作ることができます。

重要となるのは以下の区別できないという関係です。

DFA \(M(Q,\varSigma,\delta,q_0,F)\) について、状態 \(p,q \in Q\)区別できない

\(\Longleftrightarrow\) 全ての文字列\(s \in \varSigma^\ast\)について\(\delta^\ast (p,s) \in F \Leftrightarrow \delta^\ast (p,s) \in F\)

DFAの2つの状態が区別できないとことを言い換えると、2つの状態がDFAの中で同じ役割をしていると言う事ができます。DFAがどちらの状態にあってもそれから読む文字列に対して全く同じ反応を示すからです。「区別できない」という名前もここからきたのでしょう。

あるいは、区別できない2つの状態はDFAの中で無駄な状態だともいえます。なぜなら、区別できない状態の遷移を上手くいじれば2つの状態を1つにまとめてしまえることを示せる定理があるからです。いわゆるMyhill-Nerodeの定理(の一部)です。

DFA \(M(Q,\varSigma,\delta,q_0,F)\)が最小DFAである。

\(\Longleftrightarrow\) 任意の異なる状態の組\((p,q) \in Q \times Q (p \neq q)\)が区別できる。

すごくざっくり言うと、「最小DFAには"無駄な"状態が存在しない / ”無駄な”状態が無いDFAは最小」ということになります。

この定理の証明は区別できないという関係の同値類を使って遷移の具体的な張り方を示すことで行います。ここでは省略します。

というわけでDFAを最小化するには、状態のうち区別できるペアを見つけ、そのペアを一つの状態として良い感じに遷移を張ればよいということになります。以下の性質を使います。どちらも定義から明らかです。

  • \(p \in F\)\(q \notin F\)について、\(p\)\(q\)は区別できる。
  • \((p, q) \in Q \times Q\)について、ある\(c \in \varSigma\)について\(\delta(p, c)\)\(\delta(q, c)\)が区別できるならば、\(p\)\(q\)も区別できる。

これをプログラムにすると以下のようになります。

class DeterministicFiniteAutomata(Automata):
    """ 決定性有限オートマトン """
    ...
    (略)
    ...
    def minimized(self):
        """ 最小化されたDFAを返す """

        # markedとunmarkedを初期化
        # markedは区別できるペアの集合
        # unmarkedは区別できないペアの集合
        marked = set()
        unmarked = set()
        checked = set()
        for p in self.states:
            for q in self.states:
                if frozenset({p, q}) in checked or p == q:
                    continue
                if (p in self.final_states and q not in self.final_states) or \
                        (q in self.final_states and p not in self.final_states):
                    marked.add(frozenset({p, q}))
                else:
                    unmarked.add(frozenset({p, q}))
                checked.add(frozenset({p, q}))
                checked.add(frozenset({q, p}))

        # markedを限界まで増やす
        flag = True
        while flag:
            flag = False
            for p, q in unmarked:
                for s in self.alphabets:
                    # (p,q)をsで遷移させてmarkedにはいるなら(p,q)もmarkedに入る
                    if frozenset({self.transitions[p][s], self.transitions[q][s]}) in marked:
                        flag = True
                        marked.add(frozenset({p, q}))
                        unmarked.remove(frozenset({p, q}))
                        break
                if flag:
                    break

        # まず最小DFAの状態がどうなるか計算する
        states_dict = {}
        for p in self.states:
            states_dict[p] = {p}
        for p in self.states:
            for q in self.states:
                if frozenset({p, q}) in unmarked:
                    states_dict[p].add(q)
                    states_dict[q].add(p)

        # 次にstates_dictからtransitions, init_state, final_statesを作る
        states = set()
        init_state = None
        final_states = set()
        transitions = {}
        for p in self.states:
            state = frozenset(states_dict[p])
            if state in states:
                continue
            states.add(state)
            transitions[state] = {}
            for s in self.alphabets:
                # next(iter(X))は「Xから適当に一つ取る」という意味
                transitions[state][s] = states_dict[self.transitions[next(iter(state))][s]]
            if self.init_state in state:
                init_state = state
            if len(self.final_states.intersection(state)) != 0:
                final_states.add(state)
        alphabets = self.alphabets
        return DeterministicFiniteAutomata(states, alphabets, transitions, init_state, final_states)

以下のように使います。たしかに状態数が減っているのがわかります。

>>> from Automata import DeterministicFiniteAutomata as DFA
>>> a, b, c, d, e, f, g, h, = range(8)
>>> states = {a, b, c, d, e, f, g, h}
>>> alphabets = {"0", "1"}
>>> transitions = {
...     a: {"0": b, "1": f},
...     b: {"0": g, "1": c},
...     c: {"0": a, "1": c},
...     d: {"0": c, "1": g},
...     e: {"0": h, "1": f},
...     f: {"0": c, "1": g},
...     g: {"0": g, "1": e},
...     h: {"0": g, "1": c}
... }
>>> init_state = a
>>> final_states = {c}
>>> automata = DFA(states, alphabets, transitions, init_state, final_states)
>>> print(automata)

Automata
    states      : {0, 1, 2, 3, 4, 5, 6, 7}
    alphabets   : {'1', '0'}
    transitions : 0 : {'1': 5, '0': 1}
                  1 : {'1': 2, '0': 6}
                  2 : {'1': 2, '0': 0}
                  3 : {'1': 6, '0': 2}
                  4 : {'1': 5, '0': 7}
                  5 : {'1': 6, '0': 2}
                  6 : {'1': 4, '0': 6}
                  7 : {'1': 2, '0': 6}
    init_state  : 0
    final_states: {2}

>>> minimized = automata.minimized()
>>> print(minimized)

Automata
    states      : {{2}, {6}, {0, 4}, {1, 7}, {3, 5}}
    alphabets   : {'1', '0'}
    transitions : {6} : {'1': {0, 4}, '0': {6}}
                  {2} : {'1': {2}, '0': {0, 4}}
                  {0, 4} : {'1': {3, 5}, '0': {1, 7}}
                  {3, 5} : {'1': {6}, '0': {2}}
                  {1, 7} : {'1': {2}, '0': {6}}
    init_state  : {0, 4}
    final_states: {{2}}

感想

正規表現やりたいんだけどできるかな...

pythonでスクリーンショットを正しくとる(Windows)

問題

pythonの画像処理モジュールPillowにはスクリーンショットをとる関数PIL.Imagegrab.grab()があります。しかしこれはWindowsの一部の環境だとうまく動きません。具体的には、スクリーンショットをとっても左上の部分しか記録されません。

from PIL import ImageGrab
img = ImageGrab.grab()
img.save("screenshot.png")

上のコードを実行すると以下の画像が保存されます。 f:id:e-yuki67:20170212152645p:plain

この画像は本来保存されるべき画像とは異なっています。 f:id:e-yuki67:20170212152635p:plain

少し調べたところ、Windowsのスケーリング機能(デスクトップで右クリック->ディスプレイ設定で見れる)を使っているのが原因なことが分かりました。この機能は解像度が高いディスプレイでウィンドウや文字を大きく表示するために使われる機能なのですが、OSの内部ではディスプレイの解像度を小さくして座標の計算などが行われているらしく、この機能を使っているとディスプレイの解像度を取得する関数が正しい解像度を返してくれないようです。

その結果として、例えばスケーリング倍率を1.5倍にしていると解像度が1/1.5=2/3になったものとして扱われ、左上の2/3の領域しかキャプチャされません。

解決法

さらに調べると同じ問題とその解決法をPillowのgithubリポジトリissueで見つけました。そこには正しくスクリーンショットを取るコードも書いてあり、自分の環境でもうまくいきました。

issueに書いてあるコードはグローバル変数を使っているのであまりきれいなコードとは言えなかったのですが、自分で書いて自分で使うコードに使うだけだったのでそのまま使いました。ディスプレイの本当の解像度を取得する方法があるはずですが、見つけることができませんでした。

動いたコードを書いておきます。実行にはpywin32というモジュールが必要ですが、このモジュールはpipからインストールできないようなので注意が必要です。ググればインストーラーが出てきます。

import win32gui
import win32ui
import win32con
from PIL import Image

SCREEN_WIDTH = 2160
SCREEN_HEIGHT = 1440
SCREEN_SCALING_FACTOR = 1.5

def screenshot():
    """ スクリーンショット撮ってそれを(Pillow.Imageで)返す """
    window = win32gui.GetDesktopWindow()
    window_dc = win32ui.CreateDCFromHandle(win32gui.GetWindowDC(window))
    compatible_dc = window_dc.CreateCompatibleDC()
    width = SCREEN_WIDTH
    height = SCREEN_HEIGHT
    bmp = win32ui.CreateBitmap()
    bmp.CreateCompatibleBitmap(window_dc, width, height)
    compatible_dc.SelectObject(bmp)
    compatible_dc.BitBlt((0, 0), (width, height), window_dc, (0, 0), win32con.SRCCOPY)
    img = Image.frombuffer('RGB', (width, height), bmp.GetBitmapBits(True), 'raw', 'BGRX', 0, 1)
    return img

成果物

この関数を使って出来上がったのがこれです。同じサイズのスクリーンショットを複数とるのに適しています。

github.com

オートマトンをpythonで書く

(非)決定性有限オートマトンチューリングマシンより真に弱い計算能力を持つというのはどのオートマトンの教科書にも書いてあることです。このことが意味するのは、オートマトンでできることは全てチューリングマシンでできるということです。ということでチューリングマシンである普通のコンピューターでオートマトンを動かすプログラムをpythonで書きました。3か月間オートマトンについての講義を受け、その期末試験の直後にプログラムを書き始めたので、驚くほど簡単に書けました。

github.com

続きを読む