genetic algorithm to find local maxima of a function

genetic.clj

(ns genetic
  "Implementation of the algorithm described in Ch. 16, 'The New
Turing Omnibus'")

(defprotocol Chromosome
  (genes [this]
    "This chromosome's genes.")
  (value [this]
    "The value that this chromosome represents.")
  (offspring [this partner]
    "The sequence of chromosomes resulting from crossing this with
     another.")
  (mutate [this]
    "Randomly mutate."))

(defn powers-of [n]
  (iterate (partial * n) 1))

(defn transpose [& colls]
  (apply map vector colls))

(deftype Maximizer [g]
  Chromosome
  (genes [this]
    (vec g))
  (value [this]
    (let [powers (powers-of 2)]
      (/ (reduce (fn [xs [bit power]] (+ xs (* bit power)))
                 0
                 (filter (comp (complement zero?) first)
                         (transpose (reverse (genes this)) powers)))
         (nth powers (count (genes this))))))
  (offspring [this partner]
    (let [[mine theirs] (map genes [this partner])
          half-length (/ (count mine) 2)
          [a1 a2 b1 b2] (mapcat #(split-at half-length %) [mine theirs])]
      [(Maximizer. (concat a1 b2)) (Maximizer. (concat a2 b1))]))
  (mutate [this]
    (let [mine (genes this)
          flip #({0 1} % 0)]
      (Maximizer. (update-in mine [(rand-int (count mine))] flip))))

  Object
  (toString [this]
    (str (value this) "," (genes this))))

(defn rand-maximizer []
  (Maximizer. (repeatedly 6 #(rand-int 2))))

(defn rand-maximizers [n]
  (repeatedly n rand-maximizer))

(defmacro locals
  "Returns a map of locals, as keywords, to their values."
  []
  (let [locals (keys &env)]
    `(zipmap ~@(map vec [(map keyword locals) locals]))))

;;; selects the top 10 solutions (chromosomes)
;;; from those, selects the top 6 and breeds them in 3 pairs to obtain 3 new progeny
;;; mutates one of the top six solutions and one of the bottom four
;;; returns 15 solutions
(defn evolve [mystery-fn chromosomes]
  (let [fittest (->> chromosomes
                     (sort-by (comp mystery-fn value))
                     reverse
                     (take 10))
        [top bottom] (split-at 6 fittest)
        progeny (mapcat (partial apply offspring) (partition 2 top))
        mutations (map (comp mutate rand-nth) [top bottom])]
    ;; (pprint (locals))
    (vec (concat top
                 bottom
                 (take 3 progeny)
                 mutations))))

(defn evolutions [mystery-fn]
  (iterate #(evolve mystery-fn %) (rand-maximizers 15)))

(comment

  (defn f [x] (Math/sin (inc x)))

  ;;  maxima plot: plot2d(sin(x+1),[x,0,1],[gnuplot_term,dumb]);
  ;;
  ;;  sin(x+1)
  ;;       1 ++-----------+-------------+----$$$$$$$$$$$$$---------+-----------++
  ;;         +            +             $$$$$        +    $$$$$$   +            +
  ;;    0.98 ++                      $$$$                      $$$$            ++
  ;;         |                    $$$                              $$           |
  ;;    0.96 ++                $$$                                   $$$       ++
  ;;         |               $$$                                        $$      |
  ;;    0.94 ++            $$$                                            $$   ++
  ;;         |            $$                                                $$  |
  ;;    0.92 ++         $$                                                    $$+
  ;;         |        $$                                                       $$
  ;;     0.9 ++      $$                                                        ++
  ;;         |     $$                                                           |
  ;;    0.88 ++   $$                                                           ++
  ;;         |  $$                                                              |
  ;;    0.86 ++$$                                                              ++
  ;;         +$           +             +            +             +            +
  ;;    0.84 $+-----------+-------------+------------+-------------+-----------++
  ;;         0           0.2           0.4          0.6           0.8           1
  ;;                                         x

  ;; guess max(f(x)) for x in interval [0,1]
  (doseq [guess (take 30 (evolutions f))]
    (prn (-> guess first value float)))
  ;; ... 0.578125

  )