トップ > アルゴリズム > 遺伝的アルゴリズムをRで書いてみた

遺伝的アルゴリズムをRで書いてみた

アルゴリズム R Ruby

Rejectセキュリティ&プログラミングキャンプでは極度の体調不良*1により、あんまりコーディングしてる時間が取れなかったのですが、それでもなんとか自分がやろうと思っていたところができました。集合知プログラミングの11章「進化する知性」ということで、遺伝的アルゴリズムをやってみました。最初はRubyで書いてたんですが、Rubyのことをあんまり分かってないらしく、バグってるっぽい流れ*2。ということで手がちょっとは動くRで書いてみました。クラスとかを作らないと厳しい感じだったので、初めてのS4とかに挑戦してみたんですが、あんまりいいコードではない気がする。

この章でやっていることは、xとyに0から40までの一様乱数を用意して、x^2 + 3x + 2y + 5の結果を予測する関数を生成するプログラムを書こう、ということです。

関数を生成するプログラム、ということで、「解析木」というものが必要になるそうです。解析木というのはこういうようなもの。

これを下の木から評価していきます。で、どんどん木の上に渡していって、最後に一つの値を返します。で、最初はランダムな木を生成して、それを突然変異とか、交叉とかを繰り返していきます。

突然変異

突然変異は一つの木を受け取って、それをちょっと変更するという操作をするものです。サブツリーがごっそリ変わる時もあるし、全く変更しない時もあります。が、突然変異が起こる割合はあんまり高すぎないほうがいいそうです。

ただ、突然変異が起こらなすぎても、多様性が保てなくなるので、そこら辺の調整が必要なようです。

交叉

交叉というのは、二つの木を受け取って、それらを組み合わせて、新しい木を作るという操作です。あとで出てきますが、この二つの木はなるべく成功しているものが組み合わせられます。成功しているもの同士を組み合わせると、もっといいものができるだろう、という予想に基づいているようです。進化論的とかそういう感じなんでしょうか、よく分からない。

多様性

例えば、下のような関数の最大化を考えましょう。右の赤が最もよいところです。

x軸とy軸の二次元平面であると考えます。x軸の値をちょっとづつずらしていって、最大になる点というのを見ていくんですが、左の赤の最大と見なしてしまうことがあります。これを「局所最大に逹する」といいます。わりといいんだけど、一番よいというわけではない、という感じです。低次元なら、こういう図が書けるので、いいんですが、高次元だと図とか書けないので、色々な最適化の手法が考えられています。

遺伝的アルゴリズムで局所最適に陥いらないようにするための方法が、突然変異です。なんだろう、さっきの図で行くと、あるx_0の周辺でよい点を探していたんだけど、突然そこから離れた点もよいものにならないか見にいくというような感じです。交叉は近くでいい点を探しに行くような操作と見ることもできますかね。

数学的関数を再構築する

一様乱数X、Yがあった時に、x^2 + 3x + 2y + 5に適用したものと、遺伝的アルゴリズムで構築された解析木で計算した結果の絶対値の和を考えます。解析木で構築されたものが、元の関数をうまく表現できれいれば、この絶対値の和は0になります。

0になるか、規定のループ回数を終えても0にならない時終了となって、最後に木の構造が出力されます。例えばこんなのが出力される。

add
 add
  add
   if
    subtract
     p1
     add
      p1
      p0
    p1
    p1
   5
  add
   p1
   subtract
    p0
    subtract
     subtract
      subtract
       isgreater
        p0
        p0
       p0
      0
     p0
 multiply
  p0
  p0

この木をくずしていくとこんな数式が得られます。あれ、微妙に違うんだけど、絶対値の和は0になっているらしい。。。

x^2 + 4x + y + 5

元の関数がx^2 + 3x + 2y + 5なので、まあ似ているかなという感じです。どんぴしゃじゃない時もあるのかな*3。

感想

遺伝的アルゴリズムを使って、なんとかフィルタリングみたいなのをやりたいなーと思っていたんですが、それにしては時間がかかりすぎます。二重のforの中で再帰を呼びまくっているので、まあ仕方ないと言えば仕方ないのか*4。半日とか普通にかかってしまう計算量です。

なので、遺伝的アルゴリズムを使うコストに見合うものに適用するとかがいいのかなーと思います。リアルタイムでフィルタリングとかは明かに間に合いません。そういう場合は単純な計算な割りに結構精度が出るナイーブベイズとかのほうがいいようです。

  • 遺伝的アルゴリズムは(この場合は)問題の性質とかを見ずに色々な問題に対応できる
    • でも、時間はかかる
  • その他の手法は問題の性質に着目することで計算を早くしている
    • しかし、適用できるリージョンは狭くなってしまう

ということで、全ての側面でよいアルゴリズムとかはないんだよねーという当たり前の性質を再確認するということになりました。

集合知プログラミング

集合知プログラミング

以下ソースを載せていきます。

ソースコード

最新のコードはgithubに上げてあります。なんかファイルがごちゃごちゃしているけど、気にしないでくださいwww。

Rのコード

S4のところがよく分からないのですが、特にこれがないとevaluate.node関数とかが定義できないのがよく分からない。Rの基礎とプログラミング技法とかを見るといいのかなー?Rプログラミングマニュアル (新・数理工学ライブラリ 情報工学)もよい本だけど、この辺はあんまり載っていない。Rjpwikiの例もあんまりよく分からないし。。。

setGeneric("evaluate", function(this,inp) standardGeneric("evaluate"))
setGeneric("display", function(this,indent=0) standardGeneric("display"))

Rだと

  • 関数オブジェクトを自然に扱えるのが非常によい
  • データフレームや行列に任意のクラスのインスタンスをぶちこめないのが、ちょっとだめかな?
    • おかげでlistとして扱わないといけないので、コードが見づらくなってしまう
    • 評価値ともとのnodeを入れておきたいんだけど、それができない。RubyとかPythonの配列ではもちろんできるんだけど
    • よい解決方法を知っている人は教えてください

という感じですね。一長一短。ということで、Rのソースコード。

#関数のラッパークラス
setClass("fwrapper", representation(f="function", childcount="numeric",name="character"))

#何かよく分からないけど、必要っぽい
#ToDo:後で調べる
setGeneric("evaluate", function(this,inp) standardGeneric("evaluate"))
setGeneric("display", function(this,indent=0) standardGeneric("display"))

#nodeクラスの定義
setClass("node", representation(fw="fwrapper", children = "list"))

setGeneric("evaluate.node", function(this,inp) standardGeneric("evaluate.node"))
setMethod("evaluate","node",function(this,inp){
  result <- unlist(Map(function(x){evaluate(x,inp)},this@children),recursive=FALSE)
  return(this@fw@f(result))
})

setGeneric("display.node", function(this,indent) standardGeneric("display.node"))
setMethod("display","node",function(this,indent){
  cat(paste(Reduce(function(init,x){paste(init,x,sep="")},rep(" ",indent),""),this@fw@name,sep=""),fill=TRUE)
  for(c in this@children){
    display(c,indent+1)
  }
})

#paramnodeクラスの定義
setClass("paramnode", representation(idx = "numeric"))
setGeneric("evaluate.paramnode",function(this,inp){
  standardGeneric("evaluate.paramnode")
})
setMethod("evaluate","paramnode",function(this,inp){
  #PythonとRのoriginの違いを吸収するために1を足す
  return(inp[this@idx +1])
})

setGeneric("display.paramnode",function(this,indent){
  standardGeneric("display.paramnode")
})
setMethod("display","paramnode",function(this,indent){
  cat(sprintf("%sp%d",
            Reduce(function(init,x){paste(init,x,sep="")},rep(" ",indent),""),
            this@idx)
      ,fill=TRUE)
})

#constnodeクラスの定義
setClass("constnode", representation(v = "numeric"))
setGeneric("evaluate.constnode",function(this,inp){
  standardGeneric("evaluate.constnode")
})
setMethod("evaluate",signature(this="constnode"),function(this,inp){
  return(this@v)
})
setGeneric("display.constnode",function(this,indent){
  standardGeneric("display.constnode")
})
sprintf("%d%s",0,"hoge")

setMethod("display","constnode",function(this,indent){
  cat(sprintf("%s%d",
              Reduce(function(init,x){paste(init,x,sep="")},rep(" ",indent),""),
              this@v),
      fill=TRUE)
})

#関数のラッパークラスのインスタンス群
addw <- new("fwrapper",f=function(x){return(x[1]+x[2])},childcount=2,name="add")
subw <- new("fwrapper",f=function(x){return(x[1]-x[2])},childcount=2,name="subtract")
mulw <- new("fwrapper",f=function(x){return(x[1]*x[2])},childcount=2,name="multiply")

ifw <- new("fwrapper",f=function(x){
  return(ifelse(x[1] > 0,x[2],x[3]))
},childcount=3,name="if")

gtw <- new("fwrapper",f=function(x){
  return(ifelse(x[1] > x[2],1,0))
},childcount=2,name="isgreater")

flist <- c(addw,mulw,subw,ifw,gtw)

exampletree <- function(){
  return(new("node",fw=ifw,
             children=c(
               new("node",fw=gtw,children=c(
                                   new("paramnode",idx=0),new("constnode",v=3)
                                   )),
               new("node",fw=addw,children=c(
                                    new("paramnode",idx=1),new("constnode",v=5)
                                    )),
               new("node",fw=subw,children=c(
                                    new("paramnode",idx=1),new("constnode",v=2)
                                    ))
               )
             ))
}

makerandomtree <- function(pc,maxdepth=4,fpr=0.5,ppr=0.6){
  if((runif(1) < fpr) && (maxdepth > 0)){
    f <- sample(flist,1)[[1]]
    children <- Map(function(x){makerandomtree(pc,maxdepth-1,fpr,ppr)},seq(f@childcount))
    return(new("node",fw=f,children=children))
  }else if(runif(1) < ppr){
    return(new("paramnode",idx=sample(0:(pc-1),1)))
  }else{
    return(new("constnode",v=sample(0:10,1)))
  }
}

hiddenfunction <- function(x,y){
  return(x^2 + 2*y + 3*x +5)
}

buidhiddenset <- function(){
  rows <- matrix(0,200,3,byrow = TRUE)
  for(i in seq(200)){
    x <- sample(0:40,1)
    y <- sample(0:40,1)
    rows[i,] <- c(x,y,hiddenfunction(x,y))
  }
  return(rows)
}

scorefunction <- function(tree,s){
  dif <- 0
  for(i in seq(nrow(s))){
    v <- evaluate(tree,c(s[i,1],s[i,2]))
    dif = dif + abs(v-s[i,3])
  }
  return(dif)
}

mutate <- function(t,pc,probchange=0.1){
  if(runif(1) < probchange){
    return(makerandomtree(2))
  }else{
    result <- t
    if(class(t)=="node"){
      result@children <- Map(function(c){mutate(c,pc,probchange)},t@children)
    }
    return(result)
  }
}

crossover <- function(t1,t2,probswap=0.7,top=1){
  if((runif(1) < probswap) && (top != 0)){
    return(t2)
  }else{
    result <- t1
    if((class(t1)=="node") && (class(t2)=="node")){
      result@children <- unlist(Map(function(c){crossover(c,sample(t2@children,1))},t1@children),recursive=FALSE)
    }
    return(result)
  }
}

evolve <- function(pc,popsize,rankfunction,maxgen=50,mutationrate=0.1,breedingrate=0.4,pexp=0.7,pnew=0.05){
  selectindex <- function(){
    #配列のindexが0にならないように調整
    return(round(log(runif(1)) / log(pexp)) + 1)
  }
  population <- Map(function(x){makerandomtree(pc)},seq(popsize))
  for(i in seq(maxgen)){
    scores <- rankfunction(population)
    cat(scores[[1]][[1]],fill=TRUE)
    if(scores[[1]][[1]] == 0){
      break
    }
    newpop <- c(scores[[1]][[2]],scores[[2]][[2]])
    while(length(newpop) < popsize){
      if(runif(1) > pnew){
        newpop[[length(newpop) + 1]] <- mutate(crossover(scores[[selectindex()]][[2]],
                                scores[[selectindex()]][[2]],
                                probswap=breedingrate),
                      pc,probchange=mutationrate)
      }else{
        newpop[[length(newpop) + 1]] <- makerandomtree(2)
      }
    }
    population <- newpop
  }
  display(scores[[1]][[2]])
  return(scores[[1]][[2]])
}

getrankfunction <- function(dataset){
  rankfunction <- function(population){
    #任意のデータ型を入れられるのは、listのみなので、仕方なくlistでscoreを生成
    s <- c()
    p <- list()
    for(i in seq(length(population))){
      s[i] <- scorefunction(population[i][[1]],dataset)
      p[i] <- population[[i]]
    }
    l <- c()
    n <- 1
    for(i in order(s)){
      l[[n]] <- c(s[i],p[[i]])
      n <- n + 1
    }
    return(l)
  }
  return(rankfunction)
}

rf <- getrankfunction(buidhiddenset())
e <- evolve(2,500,rf)

Rubyのコード

それっぽくは動きますが、PythonやRのやつに比べると収束がいまいちすぎるので、どっかばぐってると思われます。deep copyの付近が怪しいような気がします。

require 'pp'
class Fwrapper
  attr_accessor :function
  attr_accessor :childcount
  attr_accessor :name
  def initialize(function,childcount,name)
    @function = function
    @childcount = childcount
    @name = name
  end
  def deep_copy()
    return self.class.new(Proc.new{@function},
                          Marshal.load(Marshal.dump(@childcount)),
                          Marshal.load(Marshal.dump(@name)))
  end
end 

class AbstractNode
  #あったほうが安心した
  def evaluate()
  end
  def display()
  end
  def deep_copy()
  end
end

class Node < AbstractNode

  attr_accessor :function
  attr_accessor :name
  attr_accessor :children

  def initialize(fw,children)
    @function = fw.function
    @name = fw.name
    @children = children
  end

  def evaluate(inp)
    results = @children.map{|n| n.evaluate(inp)}
    return @function.call(results)
  end

  def display(indent=0)
    puts ' '*indent + @name
    @children.each{|c|
      c.display(indent+1)
    }
  end

  def deep_copy()
    #self.childrenにもFwrapperオブジェクトが入る場合があるので、再帰的にdeep copyをやっていく
    function = @function
    childcount =  Flist.map{|x|x.childcount if x.name == @name}.compact
    name = @name
    return self.class.new(Fwrapper.new(function,childcount,name),deep_copy_children)
    #    return self.class.new(Fwrapper.new(function,childcount,name).deep_copy,deep_copy_children)
  end

  def deep_copy_children()
    return @children.map{|c|
      if c.class == Node
        c.deep_copy
      else
        c.deep_copy
        #Marshal.load(Marshal.dump(c))
      end
    }
  end
end

class Paramnode < AbstractNode
  attr_accessor :idx
  def initialize(idx)
    @idx = idx
  end
  def evaluate(inp)
    return inp[@idx]
  end
  def display(indent=0)
    puts ' '*indent + @idx.to_s
  end
  def deep_copy()
    Marshal.load(Marshal.dump(self.class.new(@idx)))
  end
end

class Constnode < AbstractNode
  attr_accessor :v
  def initialize(v)
    @v = v
  end
  def evaluate(inp)
    return @v
  end
  def display(indent=0)
    puts ' '*indent + @v.to_s
  end
  def deep_copy()
    Marshal.load(Marshal.dump(self.class.new(@v)))
  end
end

module MyMethods
  Addw = Fwrapper.new(lambda{|l|l[0]+l[1]},2,'add')
  Subw = Fwrapper.new(lambda{|l|l[0]-l[1]},2,'subtract')
  Mulw = Fwrapper.new(lambda{|l|l[0]*l[1]},2,'multiply')
  def iffunc(l)
    if l[0] > 0
      return l[1]
    else
      return l[2]
    end
  end
  module_function :iffunc

  Ifw = Fwrapper.new(self.method(:iffunc),3,'if')

  def isgreater(l)
    if l[0] > l[1]
      return 1
    else
      return 0
    end
  end
  module_function :isgreater

  Gtw = Fwrapper.new(self.method(:isgreater),2,'isgreater')

  Flist = [Addw,Mulw,Ifw,Gtw,Subw]
end

def exampletree()
  include MyMethods
  return Node.new(Ifw,[
                  Node.new(Gtw,[Paramnode.new(0),Constnode.new(3)]),
                  Node.new(Addw,[Paramnode.new(1),Constnode.new(5)]),
                  Node.new(Subw,[Paramnode.new(1),Constnode.new(2)]),
  ])
end

def makerandomtree(pc,maxdepth=4,fpr=0.5,ppr=0.6)
  include MyMethods
  if rand < fpr and maxdepth > 0
    f = Flist[rand(Flist.size)]
    children = (1..f.childcount).map{makerandomtree(pc,maxdepth-1,fpr,ppr)}
    return Node.new(f,children)
  elsif rand < ppr
    return Paramnode.new(rand(pc-1))
  else
    return Constnode.new(rand(10))
  end
end

def hiddenfunction(x,y)
  return x**2 + 2*y + 3*x +5
end

def buidhiddenset()
  rows = []
  (1..200).each{|i|
    x = rand(40)
    y = rand(40)
    rows.push([x,y,hiddenfunction(x,y)])
  }
  return rows
end

def scorefunction(tree,s)
  dif = 0
  s.each{|data|
    v = tree.evaluate([data[0],data[1]])
    dif += (v - data[2]).abs
  }
  return dif
end

def mutate(t,pc,probchange=0.1)
  if rand < probchange
    return makerandomtree(pc)
  else
    result = t.deep_copy
    if t.instance_variable_defined?("@children")
      result.children = t.children.each{|c|mutate(c,pc,probchange)}
    end
    return result
  end
end

def crossover(t1,t2,probswap=0.7,top=1)
  if rand < probswap and top == 0
    return t2.deep_copy
  else
    result = t1.deep_copy
    if t1.instance_variable_defined?("@children") && t2.instance_variable_defined?("@children")
      choice = t2.children[rand(t2.children.length)]
      result.children = t1.children.map{|c|crossover(c,choice,probswap,0)}
      return result
    else
      return result
    end
  end
end

def selectindex(pexp)
  return rand((Math::log(rand) / Math::log(pexp)).to_i + 1)
end

def myselectindex(length)
  sum = 0.0
  (1..length).each{|i|sum += 1.0 / (i**2)}
  h = Hash.new
  tmp = 0
  (1..length).each{|i|
    tmp += (1.0 / (i**2)) / sum
    h[i] = tmp
  }
  i = 0
  r = rand
  h.sort.map{|x|x[1]}.each_with_index{|item,index|
    i = index
    if item > r
      break
    end
  }
  return i 
end

def evolve(pc,popsize,rankfunction,maxgen=500,mutationrate=0.2,breedingrate=0.1,pexp=0.7,pnew=0.1)
  population = (1..popsize).map{makerandomtree(pc)}
  scores = Array.new(Array.new)
  (1..maxgen).each{
    scores = rankfunction.call(population)
    if scores[0][0] == 0.0
      break
    end
    puts scores[0][0]
    newpop = [scores[0][1],scores[1][1]]
    while newpop.length < popsize
      if rand > pnew
        l = scores.length
        newpop.push mutate(crossover(scores[myselectindex(l)][1],
                                     scores[myselectindex(l)][1],
                                     probswap=breedingrate),
                                     pc,probchange=mutationrate)
      else
        newpop.push makerandomtree(pc)
      end
    end
    population = newpop
  }
  scores[0][1].display
  return scores[0][1]
end
def getrankfunction(dataset)
  return lambda{|population| population.map{|t|[scorefunction(t,dataset),t]}.sort{|x,y|x[0] <=> y[0]}}
  #配列の配列のsortだから注意が必要。最初の要素が同じだと後のでも比較するけど、後のが比較可能でないといけない
end

*1:バス酔いから腹痛に転じた

*2:deep copyとかがうまくできてない、とか、うまく収束しないとか

*3:まあ、複雑な関数だと0になる場合が珍しいので、よいのかもしれない?

*4:何世代かに渡って改善されないならば、breakするとかそれなりにあるけど