遺伝的プログラミングをなぞってみた

　このあいだ、図書館で借りて読んだマイケル・クライトン著「プレイ―獲物―」（早川文庫SF）はおもしろくて、借りた日に一気に読んで翌日返した。

プレイ―獲物〈上〉 (ハヤカワ文庫NV)

プレイ―獲物〈下〉 (ハヤカワ文庫NV)

「プレイ」では遺伝的プログラミングによってナノマシンが勝手に進化し、人間の制御を超え（意図的なものではないが結果として）反乱を起こす。危機にさらされた技術者たちはそれに立ち向かい、辛くも生き延びる。クライトンさんのいつものパターンだ。ではあるが、不自然さもところどころあるにはありながら、とてもおもしろい小説だった。
　それで、遺伝的プログラミングというものを知りたくなった。オーム社から発行されているオライリー(O'REILLY)の「集合知プログラミング」に遺伝的プログラミングの章があったので試してみた。

集合知プログラミング

　サンプルコードはpython。pythonでのコーディングなんてしたことがなかったけど、まあなんとかなるだろうと打ち込んで試してみた。
　ところどころ誤植があって、動かない。あらら。直してみた。

kimura-shinichi-no-macbook:~ kimurashinichi$ cat -n gp.py

     1	from random import random,randint,choice

     2	from copy import deepcopy

     3	from math import log

     4	

     5	class fwrapper:

     6	  def __init__(self,function,childcount,name):

     7	    self.function=function

     8	    self.childcount=childcount

     9	    self.name=name

    10	

    11	class node:

    12	  def __init__(self,fw,children):

    13	    self.function=fw.function

    14	    self.name=fw.name

    15	    self.children=children

    16	

    17	  def evaluate(self,inp):

    18	    results=[n.evaluate(inp) for n in self.children]

    19	    return self.function(results)

    20	

    21	  def display(self,indent=0):

    22	    print (' '*indent)+self.name

    23	    for c in self.children:

    24	      c.display(indent+1)

    25	

    26	class paramnode:

    27	  def __init__(self,idx):

    28	    self.idx=idx

    29	

    30	  def evaluate(self,inp):

    31	    return inp[self.idx]

    32	

    33	  def display(self,indent=0):

    34	    print '%sp%d' % (' '*indent,self.idx)

    35	

    36	class constnode:

    37	  def __init__(self,v):

    38	    self.v=v

    39	

    40	  def evaluate(self,inp):

    41	    return self.v

    42	

    43	  def display(self,indent=0):

    44	    print '%s%d' % (' '*indent,self.v)

    45	

    46	addw=fwrapper(lambda l:l[0]+l[1],2,'add')

    47	subw=fwrapper(lambda l:l[0]-l[1],2,'subtract')

    48	mulw=fwrapper(lambda l:l[0]*l[1],2,'multiply')

    49	

    50	def iffunc(l):

    51	  if l[0]>0: return l[1]

    52	  else: return l[2]

    53	ifw=fwrapper(iffunc,3,'if')

    54	

    55	def isgreater(l):

    56	  if l[0]>l[1]: return 1

    57	  else: return 0

    58	gtw=fwrapper(isgreater,2,'isgreater')

    59	

    60	flist=[addw,mulw,ifw,gtw,subw]

    61	

    62	def exampletree():

    63	  return node(ifw,[

    64	                  node(gtw,[paramnode(0),constnode(3)]),

    65			  node(addw,[paramnode(1),constnode(5)]),

    66			  node(subw,[paramnode(1),constnode(2)]),

    67			  ]

    68		     )

    69	

    70	def makerandomtree(pc,maxdepth=4,fpr=0.5,ppr=0.6):

    71	  if random()0:

    72	    f=choice(flist)

    73	    children=[makerandomtree(pc,maxdepth-1,fpr,ppr)

    74	              for i in range(f.childcount)]

    75	    return node(f,children)

    76	  elif random()pnew:

   148	        newpop.append(mutate(

   149	                      crossover(scores[selectindex()][1],

   150	                                scores[selectindex()][1],

   151	                                probswap=breedingrate),

   152	                      pc,probchange=mutationrate))

   153	      else:

   154	        newpop.append(makerandomtree(pc))

   155	    population=newpop

   156	  scores[0][1].display()

   157	  return scores[0][1]

   158	

   159	def getrankfunction(dataset):

   160	  def rf(population):

   161	    scores=[(scorefunction(t,dataset),t) for t in population]

   162	    scores.sort()

   163	    return scores

   164	  return rf

   165	

   166	rf=getrankfunction(buildhiddenset())

   167	evolve(2,500,rf,mutationrate=0.2,breedingrate=0.1,pexp=0.7,pnew=0.1)

kimura-shinichi-no-macbook:~ kimurashinichi$ 

　これで動くようになった。よしよし。

　最後に実行したところで結果がなかなか収束しない。まあ偶然に任せてる部分があるわけだからここまで本のとおりになるわけはないとじっと待ってみる。えんえん待って、なんとか47世代目で収束した。

kimura-shinichi-no-macbook:~ kimurashinichi$ python

Python 2.5.1 (r251:54863, Jan 17 2008, 19:35:17) 

[GCC 4.0.1 (Apple Inc. build 5465)] on darwin

Type "help", "copyright", "credits" or "license" for more information.

>>> import gp

p0

p0

0/500: 11806

1/500: 5078

2/500: 4904

3/500: 4285

4/500: 2639

5/500: 1000

6/500: 1000

7/500: 1000

8/500: 1000

9/500: 1000

10/500: 1000

11/500: 1000

12/500: 1000

13/500: 1000

14/500: 1000

15/500: 1000

16/500: 1000

17/500: 1000

18/500: 1000

19/500: 1000

20/500: 1000

21/500: 1000

22/500: 1000

23/500: 1000

24/500: 1000

25/500: 1000

26/500: 1000

27/500: 1000

28/500: 1000

29/500: 1000

30/500: 1000

31/500: 1000

32/500: 1000

33/500: 1000

34/500: 1000

35/500: 1000

36/500: 1000

37/500: 821

38/500: 782

39/500: 774

40/500: 774

41/500: 497

42/500: 412

43/500: 216

44/500: 216

45/500: 216

46/500: 216

47/500: 0

subtract

 add

  if

   8

   3

   7

  add

   2

   p1

 if

  8

  subtract

   subtract

    multiply

     p0

     if

      3

      subtract

       4

       7

      3

    p1

   multiply

    p0

    p0

  p0

>>> ^D

kimura-shinichi-no-macbook:~ kimurashinichi$ 

　もちろん収束した結果は正しい関数と同等の内容だ(hiddenfunction 関数で定義している x**2+2*y+3*x+5 と同じ)。でもすごい複雑な結果。念のため、式変形して同等であることを確かめてみた。

subtract

 (add (if (8) (3) (7)) (add  (2) (p1)))

 (if

  (8)

  (subtract

   (subtract (multiply (p0) (if (3) (subtract (4) (7)) (3))) (p1))

   (multiply (p0) (p0)))

  (p0))

=

subtract

 (add (5) (p1))

 (subtract

   (subtract (multiply (p0) (-3)) (p1))

   (multiply (p0) (p0)))

=

subtract

 5+p1

 (subtract

   (-3*p0 - p1)

   p0**2)

=

5+p1 - (-3*p0 - p1 - p0**2)

=

p0**2 + 3*p0 + 2*p1 + 5 

　世代交替中は刈り込みをすると多様性が減ってかえって結果が得られにくくなるということだが、いったん収束してしまったら、やはり刈り込みして使うのがよさそうだ。

　こうして試してみると、遺伝子にムダと思える冗長な部分が大量に含まれているのもなるほどさもありなんと思う。また、途中で適合度関数の値が1000のままずっと進展がないところがあるが、もちろん突然変異や交配による変異は同じ頻度で継続しているわけである。こういうところを見ると、生物進化における形態の長期的安定化（断続平衡説）というのも、なるほどそういうことかなと思う。高校くらいの生物分野の授業にこういう実習を入れてみるとよさそうだ。