Pythonで数値シミュレーション(3) - カミキリムシが薪を喰う(要改良)
さて.今回は以前取り上げた(帰省の収穫),薪の中をカミキリムシの幼虫が食べた跡を,ランダムウォークにいくつかの規則を足して再現できるかどうかを試したものです.
結論から言うとあのパターンを再現できているとは言えませんが,とりあえずできたところまでで上げてみようと思った次第です.
履歴依存性 & 重み付き遷移
今回のシミュレーションは,前回作ったライフゲームのプログラムを基にして作成しました.
シミュレーションで再現したいことは,
・図の灰色で示した格子のように,虫は通った位置の両端も,半分だけ食べるとする
・自身が1単位時間前にいた位置には戻れないとする(Uターンはできない,はず)
・虫は,図に示すように既に食べられた格子点(白色)と半分食べた格子点(灰色),食べられていない格子点(黒色)で,その格子を選ぶ確率を変化させる
これらを実現するものとして,最後に示すコードを書くことができました.
通常のランダムウォークに,履歴依存性と重みのついたモデルとなります.
一つ前の格子点に戻らないようにするために,一回前の情報を常に保持するようにし,その格子点から現在の格子点の座標へのベクトルを考えれば,残りの3方向はこれに回転行列を掛けたものとみなせるので,扱いやすくなります(52〜62行目).
70〜85行目でmove_ruleで指定した重み付けをして,前,右,左から進路を選びます.
88〜94行目でlatticeの更新と,周辺の点を1にします.
虫の位置を常に記録するために,worm_siteを利用していますが,クリック時の扱いが難しかったです.一回集合として扱って減算を行った後に再びnumpyのndarrayに戻しています(159〜163行目,167,168行目).
以下にコードの全体を載せます.
worm_walk.py
#! /usr/bin/env python
# -*- coding:utf-8 -*-
#
# written by yuzugosho, September 2014.
from Tkinter import *
import numpy as np
import random
import sys
import time
class LifeGame:
def __init__(self, L=30, rule="1 3 6/1 2", pattern=None):
self.L = L
# 重みづけ --- 全 半 無
self.move_rule = [int(i) for i in rule.split("/")[0].split()]
# 食べる幅 歩幅 ---! 今回は使わない
self.bite_rule = [int(i) for i in rule.split("/")[1].split()]
self.lattice = np.zeros([self.L+2, self.L+2], dtype=int)
if pattern:
for x,y in pattern:
self.lattice[x,y] = 2
self.lattice[0,:] = self.lattice[self.L+1,:] = -1
self.lattice[:,0] = self.lattice[:,self.L+1] = -1
where = np.where(self.lattice==2)
self.worm_site = np.array([[x,y] for x,y in zip(where[0], where[1])])
self.worm = len(self.worm_site) # 特に使わない
def progress(self, canvas_update, update):
L = self.L
Tmax = 2000
t = 0
move_rule = self.move_rule
self.loop = True
past = self.worm_site.copy()
rn = np.random.rand
for i, x in enumerate(self.worm_site):
sett = True
while sett:
p = rn()*4
if p < 1: x[0] -= 1
elif p < 2: x[1] += 1
elif p < 3: x[0] += 1
else: x[1] -= 1
if self.lattice[x[0], x[1]] != -1:
self.worm_site[i] = [x[0], x[1]]
sett = False
self.worm_vec = self.worm_site - past
forward = np.array([[1,0], [0,1]])
right = np.array([[0,1], [-1,0]])
left = np.array([[0,-1], [1,0]])
def moveto(from_site, vec):
return [from_site + np.dot(forward, vec),
from_site + np.dot(right, vec),
from_site + np.dot(left, vec)]
def search(from_site, vec):
return [from_site + np.dot(right, vec),
from_site + np.dot(left, vec)]
while self.loop:
try:
past_lattice = self.lattice.copy()
past_worm_site = self.worm_site.copy()
# ランダムウォーク
for i, x in enumerate(self.worm_site):
canvas_update(x[0],x[1],color="white")
box = []
for v in moveto(x, self.worm_vec[i]):
if self.lattice[v[0],v[1]] == 2:
for j in range(move_rule[0]):
box.append(v)
elif self.lattice[v[0],v[1]] == 1:
for j in range(move_rule[1]):
box.append(v)
elif self.lattice[v[0],v[1]] == 0:
for j in range(move_rule[2]):
box.append(v)
if box == []:
continue
new = random.choice(box)
# latticeの更新
self.lattice[new[0],new[1]] = 2
for v in search(new, new-x):
if self.lattice[v[0],v[1]] == 0:
self.lattice[v[0],v[1]] = 1
canvas_update(v[0], v[1], color="gray")
self.worm_vec[i] = np.array(new)-x
self.worm_site[i] = new
for x in self.worm_site:
canvas_update(x[0], x[1], "red")
update()
# time.sleep(0.1)
t += 1
if t > Tmax:
self.loop = False
except KeyboardInterrupt:
print "stopped."
break
class Draw_canvas:
def __init__(self, lg, L):
self.lg = lg
self.L = L
default_size = 640
self.r = int(default_size/(2*self.L))
self.fig_size = 2*self.r*self.L
self.margin = 10
self.sub = Toplevel()
self.sub.title("Worm Walk")
self.canvas = Canvas(self.sub, width=self.fig_size+2*self.margin,
height=self.fig_size+2*self.margin)
self.c = self.canvas.create_rectangle
self.update = self.canvas.update
self.rects = dict()
for y in range(1,self.L+1):
for x in range(1,self.L+1):
if self.lg.lattice[x,y] == 2:
live = 2
else:
live = 0
tag = (x,y)
self.rects[tag] = Rect(x, y, live, tag, self)
self.canvas.pack()
def canvas_update(self, x, y, color):
v = self.rects[(x,y)]
v.root.canvas.itemconfig(v.ID, fill=color)
class Rect:
def __init__(self, x, y, live, tag, root):
self.root = root
self.x = x
self.y = y
self.live = live
if self.live == 2: color = "red"
else: color = "black"
self.ID = self.root.c(2*(x-1)*self.root.r+self.root.margin,
2*(y-1)*self.root.r+self.root.margin,
2*x*self.root.r+self.root.margin,
2*y*self.root.r+self.root.margin,
outline="#202020", fill=color, tag=tag)
self.root.canvas.tag_bind(self.ID, '
', self.pressed) def pressed(self, event):
if self.live == 2:
self.live = 0
color = "black"
remove = set([(self.x, self.y)])
tmp = self.root.lg.worm_site
y = set([(tmp[i,0], tmp[i,1]) for i in range(len(tmp))])
tmp = np.array([[i,j] for i,j in list(y - remove)])
self.root.lg.worm_site = tmp
else:
self.live = 2
color = "red"
self.root.lg.worm_site = np.append(self.root.lg.worm_site,
np.array(self.x, self.y), axis=0)
self.root.lg.lattice[self.x, self.y] = self.live
self.root.canvas.itemconfig(self.ID, fill=color)
class TopWindow:
def show_window(self, title="title", *args):
self.root = Tk()
self.root.title(title)
frames = []
for i, arg in enumerate(args):
frames.append(Frame(self.root, padx=5, pady=5))
for k, v in arg:
Button(frames[i],text=k,command=v).pack(expand=YES, fill='x')
frames[i].pack(fill='x')
self.root.mainloop()
class Main:
def __init__(self):
L = 100
worm = 10
rule = "1 3 6/1 2"
self.top = TopWindow()
pattern = [(np.random.randint(1,L+1),np.random.randint(1,L+1))
for i in range(worm)]
self.lg = LifeGame(L, rule, pattern=pattern)
self.top.show_window("Worm Walk", (('set', self.init),),
(('start', self.start),
('pause', self.pause)),
(('save', self.pr),),
(('quit', self.quit),))
def init(self):
self.DrawCanvas = Draw_canvas(self.lg, self.lg.L)
def start(self):
self.lg.progress(self.DrawCanvas.canvas_update, self.DrawCanvas.update)
def pause(self):
self.lg.loop = False
def pr(self):
import tkFileDialog
import os
if self.DrawCanvas is None:
return 1
fTyp = [('eps file', '*eps'), ('all files', '*')]
filename = tkFileDialog.asksaveasfilename(filetypes=fTyp,
initialdir=os.getcwd(), initialfile='figure_1.eps')
if filename is None:
return 0
d = self.DrawCanvas.canvas.postscript(file=filename)
def quit(self):
self.pause()
sys.exit()
if __name__ == '__main__':
app = Main()
今後の課題
最初にも述べましたが,あの特有のパターンはうまく作ることができませんでした.
ルールのパラメータの部分をもっといじってやればそれらしいものができるかもしれません.
また,今回のシミュレーションでは1つ前の格子点の位置に影響されましたが,実際の幼虫のスケールを考えると,もう1つ2つ履歴を保持しているべきかも.
そして,未だに個人的な興味以外で,この問題を考える意義が思いつかない・・・.
何かの役にたつことがあるんでしょうか.
何か提案などありましたら,コメントいただけると幸いです.