genetic programming example

gp.fs

open System
open Microsoft.FSharp.Quotations
open Microsoft.FSharp.Quotations.ExprShape
open Microsoft.FSharp.Linq.QuotationEvaluation

module Util =
  let rnd = new Random()
  let (|Range|_|) min max x =
    if min <= x && x <= max then Some() else None

module List =
  open Util
  let rndPick (xs : list<'T>) : 'T =
    rnd.Next(0, xs.Length) |> List.nth xs
  let rndAppOnce (f : 'T -> 'T) (xs : list<'T>) =
    let pos = rnd.Next(xs.Length)
    List.mapi (fun i e -> if i = pos then f e else e) xs

module GP =
  open Util

  let private varX = new Var("x", typeof<float>)

  ///ランダムに式を生成する。木の深さはdepth以下に制限される。
  let rec genExpr (depth : int) : Expr<float> = 
    let genTerm (depth : int) =
      let num = rnd.NextDouble() * 10.0
      match rnd.NextDouble() with
      | Range 0.0 0.33 -> <@ num @>
      | _        -> <@ %(Expr.Cast (Expr.Var(varX))) @>
    let genNonTerm (depth : int) =
      let rhs = depth - 1 |> genExpr
      let lhs = depth - 1 |> genExpr
      match rnd.NextDouble() with
      | Range 0.0  0.33 -> <@ %lhs + %rhs @>
      | Range 0.33 0.66 -> <@ %lhs - %rhs @>
      | _         -> <@ %lhs * %rhs @>
    match rnd.NextDouble() with
    | Range 0.0 0.3 when depth > 0 -> genNonTerm depth
    | _              -> genTerm depth

  ///突然変異を行い、新たな式を返す。
  let rec mutation (expr : Expr<float>) : Expr<float> = 
    if rnd.NextDouble() < 0.3 then
      genExpr 3
    else
      match expr with
      | ShapeCombination(a, ([_; _] as xs)) ->
        xs |> List.map Expr.Cast |> List.map mutation |> function
        | [lhs'; rhs'] ->
          RebuildShapeCombination(a, [lhs'; rhs'])
          |> Expr.Cast
        | _      -> failwith "match failure"
      | _ -> Expr.Cast expr

  ///二つの式を交配し、新たな一つの式を返す。
  let rec crossover (lhs : Expr<float>) (rhs : Expr<float>) : Expr<float> =
    let rec extract (expr : Expr<float>) : Expr<float> = 
      if rnd.NextDouble() < 0.3 then
        expr
      else
        match expr with
        | ShapeCombination(_, ([_; _] as xs)) ->
          xs
          |> List.map Expr.Cast
          |> List.rndPick
          |> extract
        | _ -> expr
    if rnd.NextDouble() < 0.3 then
      extract rhs
    else
      match lhs with
      | ShapeCombination(a, ([_; _] as xs)) ->
        xs |> List.map Expr.Cast |> List.rndAppOnce (fun x -> crossover x rhs) |> function
        | [lhs'; rhs'] ->
          RebuildShapeCombination(a, xs)
          |> Expr.Cast
        | _ -> failwith "match failure"
      | _ -> crossover lhs rhs

  ///作成する関数。
  let f x = (x ** 2.0) + (2.0 * x) + 12.0

  ///関数fとの[0, 100]の範囲の二乗誤差を求める。
  let fitness (expr : Expr<float>) : float = 
    let lambda : Expr<float -> float> =
      Expr.Lambda(varX, expr)
      |> Expr.Cast
    let lambda = lambda.Compile() ()
    seq {
      for i in 0.0..10.0..100.0 do
        yield (f i - lambda i) ** 2.0
    } |> Seq.sum

  exception private Done

  let makeNextGeneration (pool : Expr<float> []) : Expr<float> [] =
    let fitnesses =
      Array.scan (fun (i, (_, acc)) e -> 
        let f = fitness e
        if f < 0.005 then raise Done
        (i + 1, (acc, acc + 100.0 / f))) (0, (0.0, 0.0)) pool
    let (_, (_, total)) = fitnesses.[Array.length fitnesses - 1]
    total / float(Array.length pool) |> printfn "100/fitness average : %A"
    [| for x = 0 to Array.length pool - 1 do
      let piv = rnd.NextDouble() * total
      let (i, _) = Array.find (fun (_, (min, max)) -> min <= piv && piv <= max) fitnesses
      let piv'= rnd.NextDouble() * total
      let (j, _) = Array.find (fun (_, (min, max)) -> min <= piv'&& piv'<= max) fitnesses
      yield crossover pool.[i-1] pool.[j-1] |> mutation |]

  let evolve () : Expr<float> =
    let mutable pool = Array.init 1000 (fun _ -> genExpr(3))
    try
      for i = 1 to 100 do
        pool <- makeNextGeneration pool
    with
      Done -> ()
    (Array.sortBy fitness pool).[0]
    
[<EntryPoint>]
let main _ =
  let result = GP.evolve ()
  printfn "%A" result
  0

部分木をランダムに選択する部分について、最初の枝の方が選ばれやすくなっているからちょっと不味いかも。