;; web references:
;; GA archives (sourcecode)
;;   http://www.aic.nrl.navy.mil/galist/src/
;; interesting example object functions:
;;   http://www.aridolan.com/ga/gaa/gaa.html

;; TODO use all bits from each random number when generating
;; chromosomes, instead of doing % 2 on each number
;; TODO make everything pure functional (elimite set!'s), except the
;; RNG, which is inherently stateful
;; TODO add a traveling salesman problem solving individual. To test,
;; use a problem instance where all cities are on a circle.

(use srfi-1)
(use gauche.time)

(define-syntax unit-test
  (syntax-rules (=>)
    ((unit-test expr)
     (unit-test expr => #t))
    ((unit-test expr => result)
     (unless (equal? expr result)
       (format #t "failed: ~a => ~a, expected ~a\n"
               'expr expr result)))))

(define-syntax time
  (syntax-rules ()
    ((_ label proc)
     (let* ((rtc (make <real-time-counter>))
            (result (with-time-counter rtc proc)))
       (print #`",|label|: ,(time-counter-value rtc)")
       result))))

(define (print-gc)
  (let ((stat (gc-stat)))
    (format #t "(heap-size ~d total-bytes ~d)\n"
            (get-keyword :total-heap-size stat)
            (get-keyword :total-bytes stat))))

;; TODO option to initialize x-* and y-* using system timer
(define (make-combined-multiple-recursive-generator)
  (let ((x-1 151) (x-2 2711) (x-3 6719)
        (y-1 7333) (y-2 2447) (y-3 433)
        (a1 0) (a2 63308) (a3 -183326)
        (b1 86098) (b2 0) (b3 -539608)
        (m1 2147483647) (m2 2145483479))
    (lambda ()
      (let ((x (modulo (+ (* a1 x-1) (* a2 x-2) (* a3 x-3)) m1))
            (y (modulo (+ (* b1 y-1) (* b2 y-2) (* b3 y-3)) m2)))
        (set! x-3 x-2)
        (set! x-2 x-1)
        (set! x-1 x)
        (set! y-3 y-2)
        (set! y-2 y-1)
        (set! y-1 y)
        (modulo (- x y) m1)))))

(define rng (make-combined-multiple-recursive-generator))

;; (rnd)     -> integer in 0..(- (expt 2 31) 1)
;; (rnd INT) -> integer in 0..INT
;; (rnd FLT) -> float in 0..FLT
(define cmr-max (- (expt 2 31) 1))
(define (rnd . arg)
  (cond ((null? arg) (rng))
        ((exact? (car arg)) (modulo (rng) (car arg)))
        (else (* (/ (rng) cmr-max) (car arg)))))

;; using mersenne-twister
(use math.mt-random)
(define mt (make <mersenne-twister>))
;;(define mt (make <mersenne-twister> :seed (sys-time)))

(define (rnd . arg)
  (cond ((null? arg) (mt-random-integer mt cmr-max))
        ((exact? (car arg)) (mt-random-integer mt (car arg)))
        (else (* (mt-random-real0 mt) (car arg)))))

;; -----------------------------------------------------------------
(define (shuffle l)
  (map cdr
       (sort (map (lambda (x) (cons (rnd) x)) l)
             (lambda (x y) (< (car x) (car y))))))

(define (schwartzian sequence computation-fn sort-fn)
  (map car
       (sort (map (lambda (x)
                    (list x (computation-fn x)))
                  sequence)
             (lambda (a b)
               (sort-fn (cadr a) (cadr b))))))

;; -----------------------------------------------------------------
(define (make-individual-1 . keys)
  (let ((number-of-genes 32)
        (chromosome #f)
        (phenotype #f)
        (fitness 0))
    (define domain-left 0.0)
    (define domain-right 2.0)
    (define (object-function x)
      "1 - (x - 1)^2 = 2x - x^2, domain 0..2, maximum (1, 1)"
      (- (* 2 x) (* x x)))
;;     (define (object-function x)
;;       "(x - 1)^2 + 1 = x^2 - 2x +2 , domain 0..2, maxima (0, 2), (2, 2)"
;;       (+ (- (* x x) (* 2 x)) 2))

    (define (make-random-chromosome)
      (list-tabulate number-of-genes (lambda (n) (rnd 2))))

    (define (chromosome->phenotype chromosome)
      (define (chromosome->phenotype-aux chromosome sum power)
        (cond ((null? chromosome) sum)
              (else
               (chromosome->phenotype-aux (cdr chromosome)
                                          (if (= 1 (car chromosome))
                                              (+ sum (expt 2 power))
                                              sum)
                                          (- power 1)))))
      (chromosome->phenotype-aux chromosome 0 (- number-of-genes 1)))

    (define phenotype-maximum (- (expt 2 number-of-genes) 1))

    (define (phenotype->domain x)
      (+ (* (/ x phenotype-maximum)
            (- domain-right domain-left))
         domain-left))

    (define (phenotype->fitness phenotype)
      (object-function (phenotype->domain phenotype)))

    (define (chromosome->fitness chromosome)
      (object-function (phenotype->domain (chromosome->phenotype chromosome))))

    (define (update)
      (set! phenotype (chromosome->phenotype chromosome))
      (set! fitness (phenotype->fitness phenotype)))

    (define (crossover other)
      ;; perform a single-point crossover between self and other
      ;; return a list of 2 new individuals
      (let* ((pivot-point (+ 1 (rnd (- number-of-genes 1))))
             (i1 (make-individual-1
                  :chromosome (append (take chromosome pivot-point)
                                      (drop (other 'chromosome) pivot-point))))
             (i2 (make-individual-1
                  :chromosome (append (take (other 'chromosome) pivot-point)
                                      (drop chromosome pivot-point)))))
        (list i1 i2)))

    (define (print)
      ;;(format #t "(~d ~a)\n" fitness chromosome))
      (format #t "(~d ~d)\n" fitness (phenotype->domain phenotype)))

    (unit-test (phenotype->domain 0) => domain-left)
    (unit-test (phenotype->domain phenotype-maximum) => domain-right)

    (unit-test (chromosome->phenotype
                (list-tabulate number-of-genes (lambda (n) 0)))
               => 0)
    (unit-test (chromosome->phenotype
                (list-tabulate number-of-genes (lambda (n) 1)))
               => (- (expt 2 number-of-genes) 1))

    (let-keywords* keys ((init-chromosome :chromosome #f))
      (if init-chromosome
          (set! chromosome init-chromosome)
          (set! chromosome (make-random-chromosome))))
    (update)

    (lambda args
      (apply
       (case (car args)
         ((print) print)
         ((crossover) crossover)
         ((chromosome) (lambda () chromosome))
         ((fitness) (lambda () fitness)))
       (cdr args)))))

;; A Trivial 10 Variables Maximization Problem:
;; Maximize: (x1*x2*x3*x4*x5)/(x6*x7*x8*x9*x10) where (x1..x10)=[1..10]
;; The maximum is obviously when the first five variables equal 10 and
;; the last five variables equal 1. The function value for this case
;; is 100000.
(define (make-individual-2 . keys)
  (let ((number-of-genes 10)
        (chromosome #f)
        (fitness 0))
    (define domain-left 1)
    (define domain-right 10)
    (define (object-function args)
      (/ (apply * (take args 5))
         (apply * (drop args 5))))

    (define (make-random-chromosome)
      (list-tabulate number-of-genes (lambda (n) (+ 1 (rnd 10)))))

    (define (update)
      (set! fitness (object-function chromosome))
      ;; mutate
      )

    (define (crossover other)
      ;; perform a single-point crossover between self and other
      ;; return a list of 2 new individuals
      ;; TODO mutation
      (let* ((pivot-point (+ 1 (rnd (- number-of-genes 1))))
             (i1 (make-individual-2
                  :chromosome
                  (append (take chromosome pivot-point)
                          (drop (other 'chromosome) pivot-point))))
             (i2 (make-individual-2
                  :chromosome
                  (append (take (other 'chromosome) pivot-point)
                          (drop chromosome pivot-point)))))
        (i1 'mutate)
        (i2 'mutate)
        (list i1 i2)))

    (define (mutate)
      (when (< (rnd 1000000) 1000)      ; 1 promille
        (let* ((position (rnd number-of-genes))
               (allele (list-tail chromosome position))
               (diff (if (= (rnd 2) 0) 1 -1))
               (new-value-unclamped (+ (car allele) diff))
               (new-value (cond
                           ((< new-value-unclamped 1) 1)
                           ((> new-value-unclamped 10) 10)
                           (else new-value-unclamped))))
          ;;(print #`"c ,chromosome pos ,position all ,allele diff ,diff nvu ,new-value-unclamped nv ,new-value")
          (set-car! allele new-value)
          ;;(print #`"new-c ,chromosome")
          )))

    (define (dump)
      (format #t "(~d ~a)\n" fitness chromosome))

    (let-keywords* keys ((init-chromosome :chromosome #f))
      (if init-chromosome
          (set! chromosome init-chromosome)
          (set! chromosome (make-random-chromosome))))
    (update)

    (lambda args
      (apply
       (case (car args)
         ((print) dump)
         ((crossover) crossover)
         ((mutate) mutate)
         ((fitness) (lambda () fitness))
         ((chromosome) (lambda () chromosome))) ; XXX yuck
       (cdr args)))))

;; (let ((i1 (make-individual-2))
;;       (i2 (make-individual-2)))
;;   (i1 'print)
;;   (i2 'print)
;;   (let ((new (i1 'crossover i2)))
;;     ((car new) 'print)
;;     ((cadr new) 'print)))

(define (make-population size make-individual-proc)
  (let ((size size)
        (pool #f))
    (define (select-random-individual)
      (list-ref pool (rnd size)))

    ;; FIXME exhibits exponential runtimes wrt. population size
    (define (tournament-selection tournament-size)
      (car (sort (list-tabulate tournament-size
                                (lambda (i)
                                  (select-random-individual)))
                 (lambda (a b)
                   (>= (a 'fitness) (b 'fitness))))))

    (define (generation n)
      (dotimes (i n)
        (let ((new-pool '()))
          (dotimes (j (/ size 2))
            (let ((i1 (tournament-selection 3))
                  (i2 (tournament-selection 3)))
              (set! new-pool (append new-pool
                                     (i1 'crossover i2)))))
          (set! pool new-pool))
        (print-statistics i))
      ;; TODO post-run dump to file: complete population, statistics &c
      )

    (define (dump)
      (for-each (lambda (i) (i 'print)) pool))

    (define (print-statistics n)
      (let ((total (reduce + 0 (map (lambda (i)
                                      (i 'fitness))
                                    pool))))
        (format #t "(~d population ~d total ~d average ~d)\n"
                n size total (/ total size))))

    (set! pool (list-tabulate
                size (lambda (i) (make-individual-proc))))

    (lambda args
      (apply
       (case (car args)
         ((generation) generation)
         ((print) dump)
         ((print-statistics) print-statistics))
       (cdr args)))))

(let1 p (make-population 200 make-individual-2)
  (p 'generation 20))

(print-gc)

;; http://en.wikipedia.org/wiki/Roulette_wheel_selection
;; http://en.wikipedia.org/wiki/Tournament_selection
;;
;; Tournament selection runs a tournament among a few individuals and
;; selects the winner (the one with the best fitness) for crossover.
;;
;; Selection pressure can be easily adjusted by changing the
;; tournament size. If the tournament size is higher, weak individuals
;; have a smaller chance to be selected.
;;
;; Tournament selection pseudo code:
;;
;;  choose k (the tournament size) individuals from the population at random
;;  choose the best individual from pool/tournament with probability p
;;  choose the second best individual with probability p*(1-p)
;;  choose the third best individual with probability p*(p*(1-p))
;;  and so on...
;;
;; Deterministic tournament selection selects the best individual in
;; any tournament which is the same as p=1. A 1-way tournament (k=1)
;; selection is equivalent to random selection. The chosen individual
;; can be removed from the population that the selection is made from
;; if desired, otherwise individuals can be selected more than once
;; for the next generation.
;;
;; Tournament selection has several benefits, it is efficient to code,
;; works on parallel architectures and allows the selection pressure
;; to be easily adjusted.
;;
;; alternative description
;; ^^^^^^^^^^^^^^^^^^^^^^^
;; - choose n individuals
;; - copy the best one the n individuals to the next generation
;; - repeat this for population-size times
;;
;; n is called the tournament size (usually 2 or 3). Does not need
;; sorting. With n = 3 this is equal to quadratic rank selection.

;; -----------------------------------------------------------------
;; pseudo code:
;; Choose initial population
;; Evaluate each individual's fitness
;; Repeat
;;        Select best-ranking individuals to reproduce
;;        Mate pairs at random
;;        Apply crossover operator
;;        Apply mutation operator
;;        Evaluate each individual's fitness
;; Until terminating condition (see below)