Genetic programming with lexicase selection for modal problems. (see http://faculty.hampshire.edu/lspector/pubs/wk09p4-spector.pdf)

lexicase.clj

;; c) Lee Spector (lspector@hampshire.edu), 2012

;; Clojure code for genetic programming with lexicase selection for modal problems.
;; See http://faculty.hampshire.edu/lspector/pubs/wk09p4-spector.pdf
;; Written to run with Clojure 1.3.

(ns lexicase.core
  (:require [clojure.zip :as zip]))

(def target-data
  (map #(vector % (+ (* % %) % 1))
       (range -15 15 1)))

#_(def target-data
  (map #(vector % (if (< % 0)
                    (+ (* % %) % 1)
                    (* % 7)))
       (range -15 15 1)))

#_(def target-data
  (map #(vector % (if (< % -5)
                    (+ (* % %) % 1)
                    (if (< % 5)
                      (* % 7)
                      (- (* 3 % %) 3))))
       (range -15 15 1)))

(def function-table (zipmap '(+ - * pd ifz ifneg ifpos)
                            '(2 2 2 2  3   3     3)))

(defn random-function 
  []
  (rand-nth (keys function-table)))

(defn random-terminal
  []
  (rand-nth (list 'x (- (rand-int 10) 5))))

(defn random-code
  [depth]
  (if (or (zero? depth)
          (zero? (rand-int 2)))
    (random-terminal)
    (let [f (random-function)]
      (cons f (repeatedly (get function-table f)
                          #(random-code (dec depth)))))))

(defn pd
  "Protected division; returns 0 if the denominator is zero."
  [num denom]
  (if (zero? denom)
    0
    (/ num denom)))

(defn ifz [a b c] (if (zero? a) b c))

(defn ifneg [a b c] (if (neg? a) b c))

(defn ifpos [a b c] (if (pos? a) b c))

(defn errors
  [individual]
  (let [value-function (eval (list 'fn '[x] individual))]
    (map (fn [[x y]] 
           (Math/abs 
             (float (- (value-function x) y))))
         target-data)))

(defn codesize [c]
  (if (seq? c)
    (count (flatten c))
    1))

(defn at-index 
  "Returns a subtree of tree indexed by point-index in a depth first traversal."
  [tree point-index]
  (let [index (mod (Math/abs point-index) (codesize tree))
        zipper (zip/seq-zip tree)]
    (loop [z zipper i index]
      (if (zero? i)
        (zip/node z)
        (if (seq? (zip/node z)) 
          (recur (zip/next (zip/next z)) (dec i))
          (recur (zip/next z) (dec i)))))))

(defn insert-at-index
  "Returns a copy of tree with the subtree formerly indexed by
point-index (in a depth-first traversal) replaced by new-subtree."
  [tree point-index new-subtree]
  (let [index (mod (Math/abs point-index) (codesize tree))
        zipper (zip/seq-zip tree)]
    (loop [z zipper i index]
      (if (zero? i)
        (zip/root (zip/replace z new-subtree))
        (if (seq? (zip/node z))
          (recur (zip/next (zip/next z)) (dec i))
          (recur (zip/next z) (dec i)))))))

(defn mutate
  [i]
  (insert-at-index i 
                   (rand-int (codesize i)) 
                   (random-code 3)))

(defn crossover
  [i j]
  (insert-at-index i 
                   (rand-int (codesize i)) 
                   (at-index j (rand-int (codesize j)))))

(def single-thread-mode (atom false))

(defn pmapall
  "Like pmap but: 1) coll should be finite, 2) the returned sequence
   will not be lazy, 3) calls to f may occur in any order, to maximize
   multicore processor utilization, and 4) takes only one coll so far."
  [f coll]
  (if @single-thread-mode
    (map f coll)
    (let [agents (map #(agent % :error-handler (fn [agnt except] (println except))) coll)]
      (dorun (map #(send % f) agents))
      (apply await agents)
      (doall (map deref agents)))))

(defn group-with-errors
  "Returns a vector of [program total-error errors] vectors."
  [programs]
  (vec (pmapall #(let [errs (errors %)] (vector % (reduce + errs) errs)) programs)))

;; standard tournaments
(defn select
  "Returns the program of the best [program error errors] vector in a tournament
   set with the specified size."
  [prog-err-groups tournament-size]
  (let [limit (count prog-err-groups)
        tournament-set (repeatedly 
                         tournament-size 
                         #(nth prog-err-groups (rand-int limit)))]
    (first (first (sort #(< (second %1) (second %2)) (vec tournament-set))))))

(defn third [v] (nth v 2))

;; global lexicase
#_(defn select
  [prog-err-groups _tournament-size_]
  (loop [survivors prog-err-groups
         cases (shuffle (range (count (third (first prog-err-groups)))))]
    (if (or (empty? cases)
            (empty? (rest survivors)))
      (first (first survivors))
      (let [min-err-for-case (apply min (map #(nth % (first cases))
                                             (map third survivors)))]
        (recur (filter #(= (nth (third %) (first cases)) min-err-for-case)
                       survivors)
               (rest cases))))))

(defn evolve
  [popsize]
  (println "Starting evolution...")
  (loop [generation 0
         population (group-with-errors (repeatedly popsize #(random-code 5)))]
    (let [sorted (sort #(< (second %1) (second %2)) population)]
      (println "Generation:" generation)
      (println "Best total error:" (second (first sorted)))
      (println "Best errors by case:" (third (first sorted)))
      (println "Best program:" (first (first sorted)))
      (println "     Median error:" (second (nth sorted 
                                                 (int (/ popsize 2)))))
      (println "     Average program size:" 
               (float (/ (reduce + (map count (map flatten (map first population))))
                         (count population))))
      (if (< (second (first sorted)) 0.1) ;; good enough to count as success
        (println "Success:" (first (first sorted)))
        (if (> generation 100)
          (println "Failure.")
          (recur 
            (inc generation)
            (group-with-errors
              (pmapall (fn [_]
                         (let [operator (rand)
                               tsize 7]
                           (cond (< operator 0.49) (mutate (select population tsize))
                                 (< operator 0.98) (crossover (select population tsize)
                                                              (select population tsize))
                                 :else (select population tsize))))
                       (range popsize)))))))))

(dotimes [_ 30] (evolve 1000))