curve fitting.lisp

this code sole curve fitting using genetic programming

;;; Copyright (c) John Koza, All rights reserved.
;;; U.S. Patent #4,935,877.  Other patents pending.

;============================================================

;;; Kernel

(defstruct individual
  program 
  (standardized-fitness 0)
  (adjusted-fitness 0)
  (normalized-fitness 0)
  (hits 0))

(defvar *number-of-fitness-cases* :unbound
  "The number of fitness cases")

(defvar *max-depth-for-new-individuals* :unbound
  "The maximum depth for individuals of the initial
   random generation")

(defvar *max-depth-for-individuals-after-crossover* :unbound
  "The maximum depth of new individuals created by crossover")

(defvar *fitness-proportionate-reproduction-fraction* :unbound
  "The fraction of the population that will experience fitness
   proportionate reproduction (with reselection)
   during each generation")

(defvar *crossover-at-any-point-fraction* :unbound
  "The fraction of the population that will experience
   crossover at any point in the tree (including terminals)
   during each generation")

(defvar *crossover-at-function-point-fraction* :unbound
  "The fraction of the population that will experience
   crossover at a function (internal) point in the tree
   during each generation.")

(defvar *max-depth-for-new-subtrees-in-mutants* :unbound
  "The maximum depth of new subtrees created by mutation")

(defvar *method-of-selection* :unbound
  "The method of selecting individuals in the population.
   Either :fitness-proportionate, :tournament or
   :fitness-proportionate-with-over-selection.")

(defvar *method-of-generation* :unbound
  "Can be any one of :grow, :full, :ramped-half-and-half")

(defvar *seed* :unbound
  "The seed for the Park-Miller congruential randomizer.")

(defvar *best-of-run-individual* :unbound
  "The best individual found during this run.")

(defvar *generation-of-best-of-run-individual* :unbound
  "The generation at which the best-of-run individual was found.")

(defun run-genetic-programming-system
           (problem-function
            seed
            maximum-generations
            size-of-population
            &rest seeded-programs)
;; Check validity of some arguments
  (assert (and (integerp maximum-generations)
               (not (minusp maximum-generations)))
          (maximum-generations)
          "Maximum-generations must be a non-negative ~
           integer, not ~S" maximum-generations)
  (assert (and (integerp size-of-population)
               (plusp size-of-population))
          (size-of-population)
          "Size-Of-Population must be a positive integer, ~
           not ~S" size-of-population)
  (assert (or (and (symbolp problem-function)
                   (fboundp problem-function))
              (functionp problem-function))
          (problem-function)
          "Problem-Function must be a function.")
  (assert (numberp seed) (seed)
          "The randomizer seed must be a number")
;; Set the global randomizer seed.
  (setf *seed* (coerce seed 'double-float))
;; Initialize best-of-run recording variables
  (setf *generation-of-best-of-run-individual* 0)
  (setf *best-of-run-individual* nil)
;; Get the six problem-specific functions needed to 
;; specify this problem as returned by a call to
;; problem-function
  (multiple-value-bind (function-set-creator
                        terminal-set-creator
                        fitness-cases-creator
                        fitness-function
                        parameter-definer
                        termination-predicate)
      (funcall problem-function)
;; Get the function set and its associated
;; argument map
    (multiple-value-bind (function-set argument-map)
        (funcall function-set-creator)
;; Set up the parameters using parameter-definer
      (funcall parameter-definer)
;; Print out parameters report
      (describe-parameters-for-run
        maximum-generations size-of-population)
;; Set up the terminal-set using terminal-set-creator
      (let ((terminal-set (funcall terminal-set-creator)))
;; Create the population
        (let ((population
                (create-population
                  size-of-population function-set argument-map
                  terminal-set seeded-programs)))
;; Define the fitness cases using the
;; fitness-cases-creator function
          (let ((fitness-cases (funcall fitness-cases-creator))
                ;; New-Programs is used in the breeding of the
                ;; new population.  Create it here to reduce
                ;; consing.
                (new-programs (make-array size-of-population)))
;; Now run the Genetic Programming Paradigm using 
;; the fitness-function and termination-predicate provided
            (execute-generations
              population new-programs fitness-cases
              maximum-generations fitness-function
              termination-predicate function-set
              argument-map terminal-set)
;; Finally print out a report
            (report-on-run)
;; Return the population and fitness cases
;; (for debugging) 
            (values population fitness-cases)))))))

(defun report-on-run ()
  "Prints out the best-of-run individual."
  (let ((*print-pretty* t))
    (format t "~5%The best-of-run individual program ~
               for this run was found on ~%generation ~D and had a ~
               standardized fitness measure ~
               of ~D and ~D hit~P.  ~%It was:~%~S"
            *generation-of-best-of-run-individual*
            (individual-standardized-fitness *best-of-run-individual*)
            (individual-hits *best-of-run-individual*)
            (individual-hits *best-of-run-individual*)
            (individual-program *best-of-run-individual*))))

(defun report-on-generation (generation-number population)
  "Prints out the best individual at the end of each generation"
  (let ((best-individual (aref population 0))
        (size-of-population (length population))
        (sum 0.0)
        (*print-pretty* t))
    ;; Add up all of the standardized fitnesses to get average
    (dotimes (index size-of-population)
      (incf sum (individual-standardized-fitness
                  (aref population index))))
    (format t "~2%Generation ~D:  Average standardized-fitness ~
               = ~S.  ~%~
               The best individual program of the population ~
               had a ~%standardized fitness measure of ~D ~
               and ~D hit~P. ~%It was: ~%~S"
            generation-number (/ sum (length population))
            (individual-standardized-fitness best-individual)
            (individual-hits best-individual)
            (individual-hits best-individual)
            (individual-program best-individual))))

(defun print-population (population)
  "Given a population, this prints it out (for debugging) "
  (let ((*print-pretty* t))
    (dotimes (index (length population))
      (let ((individual (aref population index)))
        (format t "~&~D   ~S    ~S"
                index
                (individual-standardized-fitness individual)
                (individual-program individual))))))

(defun describe-parameters-for-run
    (maximum-generations size-of-population)
  "Lists the parameter settings for this run."
  (format t "~2%Parameters used for this run.~
              ~%=============================")
  (format t "~%Maximum number of Generations:~50T~D"
          maximum-generations)
  (format t "~%Size of Population:~50T~D" size-of-population)
  (format t "~%Maximum depth of new individuals:~50T~D"
          *max-depth-for-new-individuals*)
  (format t "~%Maximum depth of new subtrees for mutants:~50T~D"
          *max-depth-for-new-subtrees-in-mutants*)
  (format t
     "~%Maximum depth of individuals after crossover:~50T~D"
     *max-depth-for-individuals-after-crossover*)
  (format t
     "~%Fitness-proportionate reproduction fraction:~50T~D"
     *fitness-proportionate-reproduction-fraction*)
  (format t "~%Crossover at any point fraction:~50T~D"
          *crossover-at-any-point-fraction*)
  (format t "~%Crossover at function points fraction:~50T~D"
          *crossover-at-function-point-fraction*)
  (format t "~%Number of fitness cases:~50T~D"
          *number-of-fitness-cases*)
  (format t "~%Selection method: ~50T~A" *method-of-selection*)
  (format t "~%Generation method: ~50T~A" *method-of-generation*)
  (format t "~%Randomizer seed: ~50T~D" *seed*))

(defvar *generation-0-uniquifier-table*
        (make-hash-table :test #'equal)
  "Used to guarantee that all generation 0 individuals
   are unique")

(defun create-population (size-of-population function-set
                          argument-map terminal-set
                          seeded-programs)
  "Creates the population.  This is an array of size
   size-of-population that is initialized to contain individual
   records.  The Program slot of each individual is initialized
   to a suitable random program except for the first N programs,
   where N = (length seeded-programs).  For these first N
   individuals the individual is initialized with the respective
   seeded program.  This is very useful in debugging."
  (let ((population (make-array size-of-population))
        (minimum-depth-of-trees 1)
        (attempts-at-this-individual 0)
        (full-cycle-p nil))
    (do ((individual-index 0))
        ((>= individual-index size-of-population))
      (when (zerop (mod individual-index
                        (max 1 (- *max-depth-for-new-individuals*
                               minimum-depth-of-trees))))
        (setf full-cycle-p (not full-cycle-p)))
      (let ((new-program
              (if (< individual-index (length seeded-programs))
                  ;; Pick a seeded individual
                  (nth individual-index seeded-programs)
                  ;; Create a new random program.
                  (create-individual-program
                    function-set argument-map terminal-set 
                    (ecase *method-of-generation*
                      ((:full :grow) *max-depth-for-new-individuals*)
                      (:ramped-half-and-half
                        (+ minimum-depth-of-trees
                           (mod individual-index
                                (- *max-depth-for-new-individuals*
                                   minimum-depth-of-trees)))))
                    t
                    (ecase *method-of-generation*
                      (:full t)
                      (:grow nil)
                      (:ramped-half-and-half
                        full-cycle-p))))))
        ;; Check if we have already created this program.
        ;; If not then store it and move on.
        ;; If we have then try again.
        (cond ((< individual-index (length seeded-programs))
               (setf (aref population individual-index)
                     (make-individual :program new-program))
               (incf individual-index))
              ((not (gethash new-program
                             *generation-0-uniquifier-table*))
               (setf (aref population individual-index)
                     (make-individual :program new-program))
               (setf (gethash new-program
                              *generation-0-uniquifier-table*)
                     t)
               (setf attempts-at-this-individual 0)
               (incf individual-index))
              ((> attempts-at-this-individual 20)
               ;; Then this depth has probably filled up, so
               ;; bump the depth counter.
               (incf minimum-depth-of-trees)
               ;; Bump the max depth too to keep in line with new minimum.
               (setf *max-depth-for-new-individuals*
                     (max *max-depth-for-new-individuals*
                          minimum-depth-of-trees)))
              (:otherwise (incf attempts-at-this-individual)))))
    ;; Flush out uniquifier table to that no pointers
    ;; are kept to generation 0 individuals.
    (clrhash *generation-0-uniquifier-table*)
    ;; Return the population that we've just created.
    population))

(defun choose-from-terminal-set (terminal-set)
  "Chooses a random terminal from the terminal set. 
   If the terminal chosen is the ephemeral
   :Floating-Point-Random-Constant,
   then a floating-point single precision random constant
   is created in the range -5.0->5.0.
   If :Integer-Random-Constant is chosen then an integer random 
   constant is generated in the range -10 to +10."
  (let ((choice (nth (random-integer (length terminal-set))
                     terminal-set)))
    (case choice
      (:floating-point-random-constant
        ;; pick a random number in the range -5.0 ---> +5.0.
        ;; Coerce it to be single precision floating-point.
        ;; Double precision is more expensive
        ;; A similar clause to this could be used to coerce it
        ;; to double prevision if you really need
        ;; double precision.
        ;; This is also the place to modify if you need a range
        ;; other than -5.0 ---> +5.0.
        (coerce (- (random-floating-point-number 10.0) 5.0)
                'single-float))
      (:integer-point-random-constant
        ;; pick a random integer in the range -10 ---> +10.
        (- (random-integer 21) 10))
      (otherwise choice))))

(defun create-individual-program
           (function-set argument-map terminal-set