技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

2次元のランダムウォーク(python)

import numpy as np
import matplotlib.pyplot as plt

time = 2000

state = np.zeros(2)
orbit = np.zeros(2)
for t in range(time):
    speed =  1
    randomvec = np.random.random(2) - 0.5
    randomvec *= speed/np.linalg.norm(randomvec)
    state += randomvec
    orbit = np.vstack([orbit, state])

orbit = orbit.T

plt.scatter(0, 0, color="r")
plt.plot(*orbit, lw=0.5)
plt.scatter(*state, color="r")
plt.axis("equal")
plt.show()

f:id:swdrsker:20170421210847p:plain

pythonを使ってFitzHugh南雲方程式のnullclineを描く

pythonでヌルクラインを描く方法
いくつか方法があるようだけど、matplotlibとscipyを組み合わせる方法が良さそうだった。

FitzHugh南雲方程式

FitzHugh南雲方程式とは、一言でいうと神経細胞ニューロン)の挙動を数理的に表した方程式のこと。

ニューロンダイナミクス電気生理学的にモデル化したものとしてはHodgkin-Huxley方程式が最も古く、今でもニューロンのモデル研究で使われることがある。
ホジキンハクスレイ方程式とは、ホジキンとハクスレーがヤリイカの巨大軸索における神経膜の電位変化を調べ、 4変数の微分方程式で書き下したもののことで、この研究で1963年のノーベル生理学賞を受賞している。

このモデルを出発点として色々なモデルの単純化が進められていて、
その代表例がこの FitzHugh-Nagumo方程式。
FitzHugh-Nagumo方程式は、4変数のホジキンハクスレー方程式の「本質」を残したまま
2変数で表現されるように単純化したモデルだってこと。
2変数で表現することで数学的にも解析しやすくなってニューロンモデルの研究が進んだという経緯がある。

具体的には以下の式で表せられる方程式だ。
\dot{u}=c(-v+u-u^3/3+I(t))
\dot{v}=u-bv+a
uは膜電位を表し、vニューロン内部の変数である。まあvの方はどうでもよくて、膜電位によって次のニューロンに信号が伝達するかどうかが決まるからuが大事になってくる。
I(t)は外部からの入力なので、時間変化のある変数になっている。

ここで、色々パラメータがあるけれど、ここを参考に決めた。
a=0.7, b=0.8, c=10

nullclineを見る

参考にしたサイトによると、Iが定常入力の時I=0.34付近でホップ分岐(微分方程式の解が大きく変化すること)が起きるという。
せっかくなので、それを確かめてみようと思う。

ちなみに、ヌルクライン(nullcline)とは、微分方程式においてどんな感じで解が動くかなーというのを見るための曲線で、これを見るとだいたい微分方程式の解の軌跡が予想できる。
下の図では赤や緑の曲線がnullclineに相当していて青のラインが解の軌跡を表している。

外部入力:I=0.33の時

f:id:swdrsker:20170420043215p:plain:w500
f:id:swdrsker:20170420043223p:plain:w500

外部入力:I=0.34の時

f:id:swdrsker:20170420043319p:plain:w500
f:id:swdrsker:20170420043335p:plain:w500

コード

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

I = 0.33 # external input
def d_u(u, v):
    c = 10
    return c * (-v + u - pow(u,3)/3 + I)

def d_v(u, v):
    a = 0.7
    b = 0.8
    return u - b * v + a

def fitzhugh(state, t):
    u, v = state
    deltau = d_u(u,v)
    deltav = d_v(u,v)
    return deltau, deltav

# the initial state
u0 = 2.0
v0 = 1.0

t = np.arange(0.0, 20, 0.01)

y0 = [u0, v0]  # the initial state vector
y = integrate.odeint(fitzhugh, y0, t)
u_vec = y[:,0] 
v_vec = y[:,1]

plt.figure()
plt.plot(t, u_vec, label="u")
plt.plot(t, v_vec, label="v")
plt.title("membrane potential")
plt.grid()
plt.legend()

plt.figure()
plt.plot(u_vec, v_vec)
plt.xlabel("u")
plt.ylabel("v")
plt.title("FitzHugh nullcline")
padding = 0.2
umax, umin = u_vec.max() + padding,  u_vec.min() - padding
vmax, vmin = v_vec.max() + padding,  v_vec.min() - padding
U, V = np.meshgrid(np.arange(umin, umax, 0.1), np.arange(umin, vmax, 0.1))
dU = d_u(U, V)
dV = d_v(U, V)
plt.quiver(U, V, dU, dV)
plt.contour(U, V, dV, levels=[0], colors="green")
plt.contour(U, V, dU, levels=[0], colors="red")
plt.xlim([umin, umax])
plt.ylim([vmin, vmax])
plt.grid()
plt.show()

参考サイト
フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式
[SciPy-user] nullclines for nonlinear ODEs

三目並べの全ての状態数とその遷移関係を数え上げるアルゴリズム

三目並べでとりうる盤面の状態数を数えてみた。
f:id:swdrsker:20170413102529p:plain

英語ではTicTacToeと言ってボードゲームAIの入門によく使われる題材だったりする。
ゲームの状態数はよく9マス×3状態で3^9=19683通りだとか最初は9通りで次の番は8通り…だから9!=362880通りだとか言われるけれどどれも厳密な答えではない。
実は厳密な答えって探してもあんまり見つからない。そこでプログラムを組んで数え上げてみた。

# coding:utf-8
import numpy as np
from copy import deepcopy

BLANK = 0 # fix
BLACK = 1
WHITE = 2
SIMBOL = [" ","@","O",]


class Sanmoku():
    def __init__(self, n=3):
        self.graph = dict()
        self.value = dict()
        for i in range(pow(3,9)):
            self.graph[i] = []
            self.value[i] = 0
        self.make_graph()

    def display(self, number):
        board = self.num2board(number).reshape([3,3])
        for i in range(3):
            for j in range(3):
                print "|"+SIMBOL[board[i,j]],
            print "|\n"

    @staticmethod
    def num2board(number):
        if number < 0 or number > pow(3,9): return None
        number = int(number)
        map1d = np.zeros(9).astype(int)
        for i in range(9):
            map1d[i] = number%3
            number = number/3
        return map1d

    @staticmethod
    def board2num(board):
        return np.dot(board, pow(3 ,np.arange(9)))

    @staticmethod
    def is_over(board):
        """
        judge whether game is over : black wins 1 white wins -1
        """
        board2d = board.reshape([3,3])
        for i in range(3):
            if board2d[0,i] == board2d[1,i] and board2d[1,i] == board2d[2,i]:
                if board2d[0,i] == BLACK: return 1
                if board2d[0,i] == WHITE: return -1
            if board2d[i,0] == board2d[i,1] and board2d[i,1] == board2d[i,2]:
                if board2d[i,0] == BLACK: return 1
                if board2d[i,0] == WHITE: return -1
        if board2d[0,0] == board2d[1,1] and board2d[1,1] == board2d[2,2]:
            if board2d[0,0] == BLACK: return 1
            if board2d[0,0] == WHITE: return -1
        if board2d[2,0] == board2d[1,1] and board2d[1,1] == board2d[0,2]:
            if board2d[2,0] == BLACK: return 1
            if board2d[2,0] == WHITE: return -1
        return 0

    def make_graph(self):
        que = [[0,WHITE]]
        visited = dict()
        for i in self.graph.keys():
            visited[i] = False
        visited[0] = True
        while len(que) > 0:
            [number,turn] = que.pop(0)
            board = self.num2board(number)
            empty_index = [i for i, x in enumerate(board) if x == BLANK]
            if turn == BLACK: turn = WHITE
            else: turn = BLACK
            for i in empty_index:
                vboard = deepcopy(board)
                vboard[i] = turn
                new_num = self.board2num(vboard)
                self.graph[number].append(new_num)
                if not visited[new_num]:
                    visited[new_num] = True
                    reward = self.is_over(vboard)
                    if reward:
                        self.value[new_num] = reward
                        continue
                    que.append([new_num , turn])

        for i in self.graph.keys():
            if not visited[i]:
                self.graph.pop(i)
                self.value.pop(i)


if __name__=="__main__":
    import time
    start = time.time()
    sm = Sanmoku()
    end = time.time()
    print("%d [s]"%(end-start))
    print(len(sm.value))

`Sanmoku.graph`は辞書型で今の状態に対して次に遷移しうる状態が入る。
`Sanmoku.value`は辞書型で今の状態は白が勝ちか(-1)黒が勝ちか(1)が入る。

0 [s]
5478

結果として何も置いてない状態も含めて5478通りあることがわかりました!(黒先番は固定)