読者です 読者をやめる 読者になる 読者になる

正規表現の実装(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で無限リスト

 

 

pythonで無限リスト

この前pythonのジェネレータを使って作った無限リストには、一度計算した値を覚えておく機能がないため実行速度がとても遅いという問題がありました。SICPを参考にして試行錯誤した結果ラムダ式を使った遅延評価を使うことでメモ化付きの無限リストを作ることができたので書きます。

参考文献

参考にしたのは前の記事と同じでstructure and interpretation of computer programsの3.5節です。メモ化関数(memo_proc)を使った高速化が参考になりました。また後半の級数についての関数は練習問題3.53-3.61を参考にしました。

遅延評価

遅延評価とは式の値を定義されたときには計算せず、必要になってから初めて計算するという計算方法です。つまりa = f(x)と変数aを定義してもそれだけではf(x)が計算されるとは限らず、別の部分(たとえば"b = a + 1")でaが参照されて初めてf(x)が計算されます。aが参照されなければ、f(x)が計算されないでプログラムが終わることもありえます。

pythonの評価順序

pythonのプログラムを全く変更しないで遅延評価を使うことはできません。なぜなら、pythonが「変数に式の値を代入するときは式を計算する」と言う決まりを持っているためです。

def f(x):
    """ 呼ばれるとf calledとprintする関数 """
    print("f called.")
    return x + 1

a = f(0) # こう書いた時点でf(x)は計算される
print("a is", a)
b = a + 10
print("b is", b)
f called.
a is 1
b is 11

"f called"が"a is 1"よりも先に出力されていることから、f(0)が実行されているのは"b = a + 10"が実行されているときではなくa = f(0)が実行されているときだとわかります。

遅延評価の実装

そんなpythonで遅延評価を行うために、遅延評価させたい式を呼ぶゼロ引数関数を変数に代入し、値が必要になったときはその関数を呼ぶという工夫をします。pythonには無名関数を作るラムダ記法があるので、以下のように書けます。

a = lambda: f(0)     # aを「f(0)を呼び出す引数を取らない関数」とする
print("a is", a)
b = a() + 10         # aを呼び出す、すなわちf(0)の値を取得する
print("b is", b)
print("a() is", a()) # もう一度aを呼び出す
a is <function <lambda> at 0x000001FBDE38FBF8>
f called.
b is 11
f called.
a() is 1

出力の一行目から、aがf(0)の値ではなくて関数であることがわかります。2, 4行目から、aの値が参照されるとそのたびにf(0)が計算されていることがわかります。

遅延評価の組み合わせ

遅延評価される値を組み合わせて別の値を遅延評価されるようにすることもできます。

a = lambda: f(0)
b = lambda: a() + 10
c = lambda: a() + b() + 10

print("definition end.")
print("c() is", c())
definition end.
f called.
f called.
c() is 22

b()を計算するときにはa()が一回呼ばれ、c()を計算するときにはa()とb()が一回ずつ呼ばれるので、c()を計算するときにf(0)が二回呼ばれます。このことは"f called "が二回呼ばれていることから確かめられます。またa,b,cを定義しただけではf(0)は呼ばれないので、出力の最初の行は"definition end."となっています。

無限リスト

この遅延評価の仕組みを使うと、値を無限に持つリストが作れます。

簡単な無限リスト

リストが値を計算するのは先頭の要素だけで、2番目以降の要素は遅延評価する、という考え方が基本となります。つまり、

some_list = (first_value, lambda: list_tail()) # list_tail()はリストを返す

と言う形を用い、list_tail()の中でsome_listを再帰的に用いればsome_listの長さが無限になるという仕組みです。一番簡単な例を示します。

ones = (1, lambda: ones)

for _ in range(5):
    print(ones)
    print(ones[0])
    ones = ones[1]() # リストを一つ進める
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1

「onesの第二項以降はonesそのものである」と定義しているので、for文の中の"ones = ones[1]()"でonesを一つ進めているにもかかわらず、どれだけ項を進めてもones[0]で取り出されるのは1です。

ユーティリティ関数

便利のためにユーティリティ関数を定義しておきます。

def car(seq):
    """ seqの初項を返す """
    return seq[0]

def cdr(seq):
    """ seqの第2項以降を返す """
    return seq[1]()

def seq_ref(seq, n=0):
    """ seqの第n項を表示する """
    for _ in range(n):
        seq = cdr(seq)
    print(car(seq))

def seq_print(seq, n=10, end=", "):
    """ seqの最初のn項を表示する """
    for i in range(n):
        print(car(seq), end=end)
        seq = cdr(seq)
    print("")
        
def add(s, t):
    """ 数列としての和s+tを返す """
    return (car(s) + car(t), lambda: add(cdr(s), cdr(t)))

無限リストの例

# 自然数
integer = (1, lambda: add(integer, ones))
seq_print(integer)
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
# フィボナッチ数列
fib = (1, lambda: (1, lambda: add(cdr(fib), fib)))
seq_print(fib)
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 

メモ化

問題点

ここまで使ってきた遅延評価の仕組みには同じ計算を不必要に何度もしてしてしまうという問題点があります。たとえば、この例ではf(0)が2回計算されます。

a = lambda: f(0)
print(a())
print(a())
f called.
1
f called.
1

ここではf(0)の計算にかかる時間が僅かなのであまり気になりませんが、再帰的な定義などによって何回も呼ばれるときは計算の負荷が無視できなくなります。前に定義したfib()で第30項を計算するのは厳しいでしょう。

メモ化の実装

そこで一回目の呼び出しのときに値を保存しておいて、2回目以降の呼び出しでは保存した値をそのまま返すようにします。

そのために、今まで考えてきた形の遅延評価のための関数を受取り、値を保存するようにした関数を返す関数memoize()をつくります。pythonではメソッドが変数を持つことができるので、memoize()の中で新しく関数を定義してその関数自身に関数の値をresultとしてもたせます。

def memoize(func):
    """ funcに値を保存する機能を追加したものを返す """
    def temp():
        if temp.result is None:  # もしresultが計算されていなかったら...
            temp.result = func() # funcの値を計算し、保存する
        return temp.result       # 保存された値を返す
    temp.result = None
    return temp # 新しく作った関数を返す

memoize()を使うと、遅延評価を何度行っても関連づいた関数が呼ばれるのは一度だけです。

a = memoize(lambda: f(0))
print(a())
print(a())
print(a())
print(a())
f called.
1
1
1
1

メモ化付きの無限リスト

無限リストにもメモ化を組み込むために、関数を新たに一つ定義しadd()を書き直します。

def cons(a, b):
    """ 初項a, 第二項以降がbであるリストを返す"""
    return (a, memoize(b))


def add(s, t):
    """ 数列としての和s+tを返す """
    return cons(car(s) + car(t), lambda: add(cdr(s), cdr(t)))

このconsとaddを使えばメモ化付きの無限リストが作れます。

メモ化によってfibも常識的な速度で計算できます。

fib = cons(1, lambda: cons(1, lambda: add(cdr(fib), fib)))
seq_ref(fib, 1000)
70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501

無限リストの例

もういくつか無限リストの例を示します。

# 2の階乗
two_powered = cons(1, lambda: add(two_powered, two_powered))
seq_print(two_powered)
1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 
def mul(s, t):
    """ 数列としての積s*tを返す """
    return cons(car(s) * car(t), lambda: mul(cdr(s), cdr(t)))

# 階乗
fact = cons(1, lambda: mul(fact, integer))
seq_print(fact)
1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 
def sum(s):
    """ sの累積和 """
    ans = cons(0, lambda: add(s, ans))
    return ans

# 自然数の累積和
seq_print(sum(integer))
0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 
def expand(num, den, radix):
    """ x_n = num/denのradix進表示の小数点n桁目 """
    ans = cons((num * radix) // den, lambda: expand((num*radix)%den, den, radix))
    return ans

print("1/7の10進表示")
seq_print(expand(1, 7, 10))
print("1/7の小数点第1000桁")
seq_ref(expand(1, 7, 10), 1000)
print("1/3の2進表示")
seq_print(expand(1, 3, 2))
1/7の10進表示
1, 4, 2, 8, 5, 7, 1, 4, 2, 8, 
1/7の小数点第1000桁
5
1/3の2進表示
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 

幾何級数を計算する

無限といえば級数です。

\(e\) の計算

手始めに\(e\)を計算します。

\[e = \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \frac{1}{4!} + ...\]

なので、\(e\)に収束する数列approximate_e()を定義できます。

def map(s, f):
    """ sの各項にfを適用する """
    return cons(f(car(s)), lambda: map(cdr(s), f))

approximate_e = sum(map(fact, lambda x: 1/x))
print("eの近似値")
seq_print(approximate_e, end="\n")
print("第300項")
seq_ref(approximate_e, 300)
eの近似値
0
1.0
2.0
2.5
2.6666666666666665
2.708333333333333
2.7166666666666663
2.7180555555555554
2.7182539682539684
2.71827876984127

第300項
2.7182818284590455

\(e^x\) の計算

さらに一般化して、\(x\)を受け取り\(e^x\)に収束する数列を返す関数eを定義できます。

def scale(s, k):
    """ sをk倍したものを返す """
    return map(s, lambda x: k * x)

def x_powers(x):
    """ 初項init, 公比xの等比級数を返す """
    return cons(1, lambda: scale(x_powers(x), x))

def series(s):
    """ sを係数に持つ幾何級数の各項を計算する関数を返す """
    return lambda x: mul(s, x_powers(x))

def series_sum(s):
    """ sを係数に持つ幾何級数の累積和を計算する関数を返す """
    return lambda x: sum(series(s)(x))

e_coeff = map(fact, lambda x: 1/x) # e^xの係数
e = series_sum(e_coeff)
seq_ref(e(1), 30)
seq_ref(e(2), 30)
seq_ref(e(3), 30)
seq_ref(e(4), 30)
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

直接的な \(e^x\) の計算

上のようにすれば\(e^x\)の計算ができるのですが、この定義の仕方だとはっきり言って普通に定義したほうが簡単にできます。つまり"e_coeff = map(fact, lambda x: 1/x)"などと定義しなくても

def fact_func(n):
        return n * fact_func(n-1) if n > 1 else 1

def e(x, n):
    ans = 0
    for i in range(n):
        ans += x ** i / fact_func(i)
    return ans

print(e(1, 30))
print(e(2, 30))
print(e(3, 30))
print(e(4, 30))
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

とできるのです。こうした方が読みやすい上に遅延評価がどうだとか考える必要もありません。このままだと遅延評価の意味が無いので、\(e^x\)をもう少し面白い方法で定義し直します。

\(e^x\) の陰な定義

注目するのは、「\(e^x\)積分すると定数部分を除いて\(e^x\)に等しい」という数学的な事実です。幾何級数に対して積分を行うと項が一つ後ろにずれるので、もと級数積分した級数の間に関係があれば、級数を陰に定義することができます。

まず、幾何級数の係数を受け取ってそれを積分した級数の係数を返す関数を書きます。

def integrate_series(s):
    """ 級数sを積分した級数を返す、ただし定数項は(不定なので)除く """
    return mul(s, map(integer, lambda x: 1/ x))

seq_print(fact)
e_coeff = map(fact, lambda x: 1/x) # e^xの係数
print("e^xの係数")
seq_print(e_coeff, 7)
print("e^xを積分した係数(初項を除く)")
seq_print(integrate_series(e_coeff), 7)
1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 
e^xの係数
1.0, 1.0, 0.5, 0.16666666666666666, 0.041666666666666664, 0.008333333333333333, 0.001388888888888889, 
e^xを積分した係数(初項を除く)
1.0, 0.5, 0.16666666666666666, 0.041666666666666664, 0.008333333333333333, 0.0013888888888888887, 0.0001984126984126984, 

\(e^x\)積分したときの係数が\(e^x\)の係数と初項を除いて同じことがわかります。

このことを用いて、e_coeffを再帰的に(陰に)定義します。

# e_coeffは初項が1で、それ以降は「e_coeffを積分した級数から初項を除いたもの」である
e_coeff = cons(1, lambda: integrate_series(e_coeff))
e = series_sum(e_coeff)
seq_ref(e(1), 30)
seq_ref(e(2), 30)
seq_ref(e(3), 30)
seq_ref(e(4), 30)
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

陰に定義できる他の関数

級数の陰な定義の便利な点は\(e^x\)以外の関数を定義するときにも使える点です。たとえばsinとcosは次のように定義できます。

from math import pi

cos_coeff = cons(1, lambda: scale(integrate_series(sin_coeff), -1))
sin_coeff = cons(0, lambda: integrate_series(cos_coeff))
cos = series_sum(cos_coeff)
sin = series_sum(sin_coeff)

print("cos")
seq_ref(cos(pi / 6), 30)
seq_ref(cos(pi / 4), 30)
seq_ref(cos(pi / 3), 30)
print("sin")
seq_ref(sin(pi / 6), 30)
seq_ref(sin(pi / 4), 30)
seq_ref(sin(pi / 3), 30)
cos
0.8660254037844386
0.7071067811865475
0.5000000000000001
sin
0.49999999999999994
0.7071067811865475
0.8660254037844386

級数の積

ところで、\({sin}^2(x) + {cos}^2(x) = 1\)\(x\)がいくらであってもいつでも成り立つはずです。このことを確認してみます。

そのためには\({sin}^2, {cos}^2\)級数展開を求める必要があるので、級数の掛け算をする関数mul_coeff()を定義します。mul_coeff()の計算には

\[(a_0+a_1x+a_2x^2+\cdots)(b_0+b_1x+b_2x^2+\cdots) = a_0b_0 + a_0(b_1x+b_2x^2+\cdots) + (a_1x+a_2x^2+\cdots)(b_0+b_1x+b_2x^2+\cdots)\]

という関係を利用します。

def mul_coeff(s, t):
    """ sとtの多項式としての積 """
    return cons(car(s) * car(t), lambda: add(scale(cdr(t), car(s)), mul_coeff(cdr(s), t)))

maybe_one_coeff = add(mul_coeff(sin_coeff, sin_coeff), mul_coeff(cos_coeff, cos_coeff)) 
maybe_one = series_sum(maybe_one_coeff)

seq_ref(maybe_one(0.1), 10)
seq_ref(maybe_one(0.5), 10)
seq_ref(maybe_one(1.0), 10)
seq_ref(maybe_one(2.0), 10)
seq_ref(maybe_one(3.0), 10)
seq_ref(maybe_one(4.0), 10)
seq_ref(maybe_one(10.0), 10)
1.0
1.0
1.0
1.0000000000000004
0.9999999999999993
0.9999999999999858
0.9999999999936162

\({sin}^2(x) + {cos}^2(x)\)がほぼ1になっていることが確認できます。

級数の割り算

掛け算ができたので割り算もやります。割り算は掛け算より複雑です。

級数の逆数

まず多項式の逆数を考えます。\(S\)級数として

\[ 1 \div S = X \]

を求めることを考えます。

\(S=s_c + S_r\)(\(s_c\)は定数、\(S_r\)級数)として

\[ S \cdot X = 1\] \[ (s_c + S_r) \cdot X = 1\] \[ s_c \cdot X + S_r \cdot X = 1\] \[ X = \frac{1}{s_c}(1 - S_r \cdot X) \]

となります。\(S_r\)には定数項が含まれないので、最後の式の右辺にある\(S_r \cdot X\)にも定数項が含まれないことがポイントです。これによって最初にXの定数項を計算し、それ以外の項は\(-S_r \cdot X\)の部分から逐次に計算することができます。

コードにします。

def reciprocal(s):
    """ 級数としての逆数1/sを返す """
    return cons(1/car(s), lambda: scale(mul_coeff(cdr(s), reciprocal(s)), -1))

あっているか確認するために\(sec(x)\)を計算してみます。

sec_coeff = reciprocal(cos_coeff)
sec = series_sum(sec_coeff)

maybe_one_coeff = mul_coeff(cos_coeff, sec_coeff) 
maybe_one = series_sum(maybe_one_coeff)

print("cosの係数")
seq_print(cos_coeff, 10)
print("secの係数")
seq_print(sec_coeff, 10)
print("maybe_oneの係数")
seq_print(maybe_one_coeff, 10)
cosの係数
1, -0.0, -0.5, 0.0, 0.041666666666666664, -0.0, -0.0013888888888888887, 0.0, 2.4801587301587298e-05, -0.0, 
secの係数
1.0, 0.0, 0.5, -0.0, 0.20833333333333334, -0.0, 0.08472222222222223, -0.0, 0.034350198412698416, -0.0, 
maybe_oneの係数
1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 

maybe_oneの計算はせずに係数を表示させました。cosとsecは係数に値を持っているにも関わらずそれらの積は係数が0なことから、たしかにsecがcosの逆数なことがわかります。

一般の級数の割り算

逆数の計算ができれば一般の割り算も簡単です。

def div_coeff(s, t):
    """ 級数としての逆数s/tを返す """
    return mul_coeff(s, reciprocal(t))

\(tan(x) = sin(x) \div cos(x)\)で確認します。

tan_coeff = div_coeff(sin_coeff, cos_coeff)
tan = series_sum(tan_coeff)

seq_ref(tan(0.0), 30)
seq_ref(tan(pi/6), 30)
seq_ref(tan(pi/4), 30)
seq_ref(tan(pi/3), 30)
0.0
0.5773502691896235
0.9999999992094688
1.7320428394905787

感想

SICPではこのあとにエイトケンのΔ2乗加速法を使って円周率を求める話があるのですがこの記事には書きませんでした。

疲れた。

ジェネレータで遊ぶ(python)

 

pythonにはジェネレータという組み込み型が存在し、値を逐次に計算する配列として使うことができます。ジェネレータを使って遊んだので書きます。

参考文献

structure and interpretation of computer programsの3.5節を参考にしました。

無限ジェネレーターに関するユーティリティ関数

この文章では無限に値を返すジェネレーターを扱うので、最初に便利な関数を定義しておきます。

# gの最初のn項を表示する
def g_print(g, n=10):
    for i in range(n):
        print(next(g))

# gをn項飛ばす
def g_delayed(g, n):
    for _ in range(n):
        next(g)
    return g

# gの第n項を表示する
def g_ref(g, n):
    return next(g_delayed(g, n))

無限ジェネレータ

pythonではジェネレータを使って無限に値を返し続けるものを書けます。

# 自然数を0から返す
def integers():
    i = 0
    while True:
        yield i
        i += 1
>>> g_print(integers())
0
1
2
3
4
5
6
7
8
9

同様に少し凝ったものも書けます。

# 5の倍数でない自然数を小さい順に返す
def no_five():
    i = 0
    while True:
        if i % 5 != 0:
            yield i
        i += 1
>>> g_print(no_five())
1
2
3
4
6
7
8
9
11
12

このno_fiveではローカル変数iを1づつ増やすことで全ての自然数をチェックしています。

no_fiveで行っている「ローカル変数に無限に自然数を代入して何かする」ということはまさにintegersがfor文で使われたときに行うことなので、次のように書くことができます。

def no_five():
    for i in integers() # while文の代わりにfor文とintegersを使った
        if i % 5 != 0:
            yield i

integersからno_fiveが作られる関係は、高階関数を使ってまとめることができます。

# ジェネレータg(無限でも良い)からcondを満たすものだけを抜き取る
def g_filter(g, cond):
    for i in g:
        if cond(i): # "i % 5 != 0" をcondに一般化した
            yield i

no_five = g_filter(integers(), lambda x: x % 5 != 0)

なおpythonの組み込み関数にfilterがありますが、無限ジェネレータに対しては上手く動きません。要素をすべて読みに行ってしまうようです。

g_filterと同様に、ジェネレータに対する高階関数を定義できます。

# gの全要素にfを適用する
def g_map(g, f):
    for i in g:
        yield f(i)

# gの要素とhの要素を足す
def g_add(g, h):
    for a, b in zip(g, h):
        yield a + b

# gの要素とhの要素をかける
def g_mul(g, h):
    for a, b in zip(g, h):
        yield a * b

# gの要素の累積和
def g_sum(g):
    sum = 0
    for i in g:
        sum += i
        yield sum

これらを使って様々なジェネレータを定義できます。

# 偶数全部
evens = g_map(integers(), lambda x: 2*x)

# 7で割ると3余る偶数全部
Nseven_plus_three = g_filter(evens, lambda x: x % 7 == 3)

# 7で割ると3余る偶数に自然数を足したもの
Nseven_plus_three_plus_N = g_add(Nseven_plus_three, integers())

再帰ジェネレータ

ジェネレータは関数の一種なので、再帰的に定義することができます。

def ones():
    yield 1
    yield from ones() # 次の要素は自分自身から取る

定義が言っているのは「onesは最初1を返し、それ以降はonesが返すものを最初から返す」ということです。つまり、1を無限に返します。

>>> g_print(ones())
1
1
1
1
1
1
1
1
1
1

ただしこの例では、onesはyieldするたびに新しくonesを作成するので、onesはメモリーを無限に食います。そのためこのonesは実用的ではありません。

もう少し面白い例がこれです。

def fib():
    yield 1
    yield 1
    yield from g_add(g_delayed(fib(), 1), fib())

fibの定義を解釈すると、

  1. fibの最初は1
  2. その次は1
  3. その次以降は「fibと、fibを一つずらしたジェネレータの和」というジェネレータの返す値を返す

3.が言っているのは\(x_n = x_{n-1} + x_{n-2}\)ということなので、fibはフィボナッチ数列を表します。

>>> g_print(fib())
1
1
2
3
5
8
13
21
34
55

メモ化していないので低速ですが、ifを使わないフィボナッチ数列の定義というのは珍しいと思います。

再帰的に定義されたジェネレータの例をいくつかあげます。他では見ないような形で再帰が使われているので面白いです。

# 1, 2, 3, 4, 5, ...
def integers():
    yield 1
    yield from g_add(ones(), integers())

# 1, 2, 4, 8, 16, ...
def two_powered():
    yield 1
    yield from g_add(s(), s())

# 1, 2, 6, 24, 120, 720, ...
def fact():
    yield 1
    yield from g_mul(g_delayed(integers(), 2), fact())

感想

本当はこのあと再帰ジェネレータをメモ化して速くして終わる予定だったんだけどメモ化ができなかったのでここで終わる。こういう細かいところはpythonじゃ無理だ。

オートマトンを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}}

感想

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

凸包と記号摂動法

 

参考文献

図書館でたまたま目に入ったというだけの理由で計算幾何学の本を読みました。この本は説明のわかりやすさもさることながら、文章がうまくて読んでいて楽しかったです。

本の最初に乗っていたアルゴリズムを実装したので解説してみます。コードはgithubにあります。

考える問題は二次元平面上の点の集合\(X\)を入力として、\(X\)の凸法を求める問題です。

逐次追加法

この問題を解く逐次追加法と呼ばれるアルゴリズムがあります。\(x\)座標が小さい順に考慮する点を増やし、各段階で\((x_1, \cdots, x_i)\)に対する凸包\(Y\)を作る、という名前の通りのアルゴリズムです。

入力: 点の集合\(X\)

出力:\(X\)の凸包\(Y\)

アルゴリズム :

  1. \(X\)\(x\)座標で昇順にソートし、ソート結果を\((x_1, \cdots, x_n)\)とする。
  2. \(x_1, x_2, x_3\)を左回りになるように並べたものを\(Y\)とする。
  3. \(i=4\)とする。
  4. \(x_p\)を、「\(x_i\)から見える\(X\)の点の中で一番右にある点」とする。
  5. \(x_q\)を、「\(x_i\)から見える\(X\)の点の中で一番左にある点」とする。
  6. \(Y\)に属する点で、\(x_p\)\(x_q\)の間にあるものを\(Y\)から削除する。
  7. \(Y\)\(x_p\)の直後に\(x_i\)を追加する。
  8. \(i = n\)なら\(Y\)を出力して終了。
  9. \(i = i + 1\)として、\(4\)に戻る。

ただし、点の集合\(X\)\(X\)の外部の点\(p\)に対して、

 \(x \in X\)\(p\)から見える\(\Leftrightarrow\)\(p\)\(x\)を結ぶ直線が\(X\)の内部を通らない

とする。

逐次追加法のステップ3, 4では、点の集合の中で一番右端(あるいは左端)にある点を調べるのに以下の関係を使います。

\[F(a,b,c)=\left| \begin{array}{ccc} 1 & a_x & a_y \\ 1 & b_x & b_y \\ 1 & c_x & c_y \end{array} \right|\] について、

\(F(a,b,c) \gt 0\)ならば3点\((a, b, c)\)がこの順番に左に折れている。

\(F(a,b,c) \lt 0\)ならば3点\((a, b, c)\)がこの順番に右に折れている。

退化

逐次追加法は\(X\)の点が常識的に位置している場合には正しい凸包を計算できますが、例外的な場合には計算に失敗することがあります。

例外的な場合とは、ステップ4,5で\(F(p,q,r)=0\)となる場合です。これは3点が一直線上にある、または3点の座標が全く同じ場合に相当します。このような、一般的な方法で解くことができない例外的な(たちの悪い)状態を「退化」といいます。

退化なんて例外的なんだからほっておいてもいいだろうと考えたくなりますが、そういうわけにも行きません。なぜなら、コンピュータで扱う数値は有限の精度しか持たない以上一致することはあり得るうえに、退化が起きた結果正しい答えとまるっきり違った解答をしてしまう可能性があるためです。もしこれがソートの問題なら、退化が起きたとしても問題にはならなかったでしょう(同じ数値の点はどう並べても良いため)。

記号摂動法

退化の問題を解決する方法の一つが記号摂動法です。記号摂動法の考え方は、座標が全く同じだとまずいのだから、すべての点を微小量ずつランダムにずらして同じ値を持つ点をなくしてしまえば良い、というものです。微小量ずらすと言っても本当に値をずらしてしまうとずらさずにうまく行っていた処理がうまくいかなくなる可能性があるので、値をずらすことはせずに、微小量は多項式の項として追加するという方法を取ります。

具体的な手順は以下のようになります。

入力された点の集合\[X=\{(x_1, y_1), (x_2, y_2), \cdots, (x_n, y_n)\}\] に対し、十分大きな整数\(M\)を用いて\(\varepsilon\)多項式\(x'_i, y'_i\)\((i = 1, \cdots, n)\)\[(x'_i, y'_i) =(x_i+\varepsilon^i, y_i + \varepsilon^{Mi})\] と定義する。そして \[X'=\{(x'_1, y'_1), (x'_2,y'_2), \cdots, (x'_n, y'_n)\}\] についてアルゴリズムを適用する。 ただし、\(\varepsilon\)多項式の符号は、多項式\(0\)で無い項のうち一番次数の小さなものの符号とする。

このようにすることで、\(F(a, b, c)\)の符号の計算について以下のような性質が成り立ちます。

  • 記号摂動を使えば退化は起きない。
    • 退化が起こるとき\(F(a, b, c)\)\(\varepsilon\)多項式として\(0\)となるが、\(a_x,a_y,b_x,b_y,c_x,c_y\)は全て違う\(\varepsilon\)の次数を持つことからそのようなことは起こり得ない。
  • 記号摂動を使わないでも退化が生じていない場合、記号摂動を用いても同じ結果が得られる。
    • 退化が生じていない場合には\(F(a, b, c)\)の符号が定数項の符号で判定されるため。
  • 記号摂動を使わないと退化が生じる場合、記号摂動によって\(F(a, b, c)\)の符号が一意に定まり、かつ幾何学的な性質を満たす。
    • 全ての\(i\)について\(x'_i\)\(y'_i\)に順序が定まっているため。

このような性質によって、退化を場合分けに頼ることなく一気に解決することができます。

おまけ: ギフトラッピング法

凸法を求めるアルゴリズムは逐次追加法の他にもギフトラッピング法と呼ばれるアルゴリズムがあります。

\(X\)の中で一番右端にある点に紐の片方の端をくくりつけ、紐のもう片方の端を持って左回りに回していったときに紐が触れる点と同じ順番で\(Y\)に点が追加されていきます。3次元で同じことをやると箱を包装紙で包むやり方に似るのでギフトラッピング法と名づけられたのでしょう。

入力: 点の集合\(X\)

出力:\(X\)の凸包\(Y\)

アルゴリズム :

  1. \(X\)の中で\(x\)座標が一番大きいものを\(p\)とする。
  2. \(p_{1}=p\)とする。
  3. \(q\)を、「\(X-\{p\}\)に属する点\(q\)のうち、\(\vec{pq}\)\(y\)軸正方向のなす角度が最小になる点」とする。
  4. \(Y = \{ \vec{pq} \}\)とする。
  5. \(U = X\)とする。
  6. \(r\)を、「\(U\)に属する点\(r\)のうち、\(\vec{pq}\)\(\vec{q r}\)のなす角度が最小になる点」とする。
  7. \(Y\)\(\vec{q r}\)を追加する。
  8. \(r\)\(p_1\)だったら、\(Y\)を出力して終了。
  9. そうでなかったら、\((p,q) = (q,r)\)として6.に戻る。

このアルゴリズムのステップ6では3点の関係を調べるのに以下の関係を使います。

平面上の\(3\)\(p, q, r\)について、\(\vec{pq}\)\(\vec{q r}\)のなす角度を\(\theta\)とすると \[cos\theta = \frac{\vec{pq} \cdot \vec{q r}}{|\vec{pq}||\vec{q r}|}\]

ステップ6では\(\vec{pq}\)が凸包の辺であることから\(0 \le \theta \le \pi\)であり、かつこの範囲において\(cos\theta\)は単調減少なことから、\(\theta\)の最小値を求めるには\(cos\theta\)の最大値を求めればいいことがわかります。

\(p=q\)または\(q=r\)のとき\(cos\theta=\infty\)となるため、このアルゴリズムでも退化が起きます。しかしこの退化は逐次追加法で起こった退化よりたちが悪いです。\(cos\theta\)は割り算を含むので、記号摂動法では計算できないためです。

この退化の解決法は本に載っておらず、また自分で思いつくこともできなかったので書けませんが、単純で簡単なアルゴリズムが常に最良であるわけではないという例だと思ったので書きました。

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

続きを読む