curve fitting.lisp

this code sole curve fitting using genetic programming

          (new-female (list (copy-tree female))))
      ;; Get the pointers to the subtrees indexed by male-point
      ;; and female-point
      (multiple-value-bind (male-subtree-pointer male-fragment)
          (get-subtree (first new-male) new-male male-point)
        (multiple-value-bind
          (female-subtree-pointer female-fragment)
            (get-subtree
              (first new-female) new-female female-point)
          ;; Modify the new individuals by smashing in the
          ;; (copied) subtree from the old individual.
          (setf (first   male-subtree-pointer) female-fragment)
          (setf (first female-subtree-pointer) male-fragment)))
      ;; Make sure that the new individuals aren't too big.
      (validate-crossover male new-male female new-female))))

(defun count-crossover-points (program)
  "Counts the number of points in the tree (program).  
   This includes functions as well as terminals."
  (if (consp program)
      (+ 1 (reduce #'+ (mapcar #'count-crossover-points
                               (rest program))))
      1))

(defun max-depth-of-tree (tree)
  "Returns the depth of the deepest branch of the
   tree (program)."
  (if (consp tree)
      (+ 1 (if (rest tree)
               (apply #'max
                      (mapcar #'max-depth-of-tree (rest tree)))
               0))
      1))

(defun get-subtree (tree pointer-to-tree index)
  "Given a tree or subtree, a pointer to that tree/subtree and
   an index return the component subtree that is numbered by
   Index.  We number left to right, depth first."
  (if (= index 0)
      (values pointer-to-tree (copy-tree tree) index)
      (if (consp tree)
          (do* ((tail (rest tree) (rest tail))
                (argument (first tail) (first tail)))
               ((not tail) (values nil nil index))
            (multiple-value-bind
                (new-pointer new-tree new-index)
                (get-subtree argument tail (- index 1))
              (if (= new-index 0)
                  (return
                    (values new-pointer new-tree new-index))
                  (setf index new-index))))
          (values nil nil index))))

(defun validate-crossover (male new-male female new-female)
  "Given the old and new males and females from a crossover
   operation check to see whether we have exceeded the maximum
   allowed depth.  If either of the new individuals has exceeded
   the maxdepth then the old individual is used."
  (let ((male-depth   (max-depth-of-tree (first new-male)))
        (female-depth (max-depth-of-tree (first new-female))))
    (values
      (if (or (= 1 male-depth)
              (> male-depth
                 *max-depth-for-individuals-after-crossover*))
          male
          (first new-male))
      (if (or (= 1 female-depth)
              (> female-depth
                 *max-depth-for-individuals-after-crossover*))
          female
          (first new-female)))))

(defun crossover-at-function-points (male female)
  "Performs crossover on the two programs at a function
   (internal) point in the trees."
  ;; Pick the function (internal) points in the respective trees 
  ;; on which to perform the crossover.
  (let ((male-point
          (random-integer (count-function-points male)))
        (female-point
          (random-integer (count-function-points female))))
    ;; Copy the trees because we destructively modify the new
    ;; individuals to do the crossover and Reselection is
    ;; allowed in the original population. Not copying would
    ;; cause the individuals in the old population to
    ;; be modified.
    (let ((new-male   (list (copy-tree male)))
          (new-female (list (copy-tree female))))
      ;; Get the pointers to the subtrees indexed by male-point
      ;; and female-point
      (multiple-value-bind (male-subtree-pointer male-fragment)
          (get-function-subtree
            (first new-male) new-male male-point)
        (multiple-value-bind
          (female-subtree-pointer female-fragment)
            (get-function-subtree
              (first new-female) new-female female-point)
          ;; Modify the new individuals by smashing in
          ;; the (copied) subtree from the old individual.
          (setf (first   male-subtree-pointer) female-fragment)
          (setf (first female-subtree-pointer) male-fragment)))
      ;; Make sure that the new individuals aren't too big.
      (validate-crossover male new-male female new-female))))

(defun count-function-points (program)
  "Counts the number of function (internal) points
   in the program."
  (if (consp program)
      (+ 1 (reduce #'+ (mapcar #'count-function-points
                               (rest program))))
      0))

(defun get-function-subtree (tree pointer-to-tree index)
  "Given a tree or subtree, a pointer to that tree/subtree and
   an index return the component subtree that is labeled with
   an internal point that is numbered by Index.  We number left
   to right, depth first."
  (if (= index 0)
      (values pointer-to-tree (copy-tree tree) index)
      (if (consp tree)
          (do* ((tail (rest tree) (rest tail))
                (argument (first tail) (first tail)))
               ((not tail) (values nil nil index))
            (multiple-value-bind
              (new-pointer new-tree new-index)
                (if (consp argument)
                    (get-function-subtree
                      argument tail (- index 1))
                    (values nil nil index))
              (if (= new-index 0)
                  (return
                    (values new-pointer new-tree new-index))
                  (setf index new-index))))
          (values nil nil index))))

(defun mutate (program function-set argument-map terminal-set)
  "Mutates the argument program by picking a random point in
   the tree and substituting in a brand new subtree created in
   the same way that we create the initial random population."
  ;; Pick the mutation point.
  (let ((mutation-point
          (random-integer (count-crossover-points program)))
        ;; Create a brand new subtree.
        (new-subtree
          (create-individual-program
            function-set argument-map terminal-set
            *max-depth-for-new-subtrees-in-mutants* t nil)))
    (let ((new-program (list (copy-tree program))))
      (multiple-value-bind (subtree-pointer fragment)
          ;; Get the pointer to the mutation point.
          (get-subtree (first new-program)
                       new-program mutation-point)
        ;; Not interested in what we're snipping out.
        (declare (ignore fragment))
        ;; Smash in the new subtree.
        (setf (first subtree-pointer) new-subtree))
      (values (first new-program) new-subtree))))

(defun park-miller-randomizer ()
  "The Park-Miller multiplicative congruential randomizer
   (CACM, October 88, Page 1195).  Creates pseudo random floating
   point numbers in the range 0.0 < x <= 1.0.  The seed value
   for this randomizer is called *seed*, so you should
   record/set this if you want to make your runs reproducible."
  #+Lucid (unless (typep *seed* 'integer) (setq *seed* (round *seed*)))
  (assert (not (zerop *seed*)) () "*seed* cannot be zero.")
  (let ((multiplier #+Lucid 16807 #-Lucid 16807.0d0);16807 is (expt 7 5)
        (modulus #+Lucid 2147483647 #-Lucid 2147483647.0d0))
             ;2147483647 is (- (expt 2 31) 1)
    (let ((temp (* multiplier *seed*)))
      (setf *seed* (mod temp modulus))
      ;;Produces floating-point number in the range
      ;;  0.0 < x <= 1.0
      (#+lucid float #-lucid progn (/ *seed* modulus)))))

(defun random-floating-point-number (n)
  "Returns a pseudo random floating-point number
                in range 0.0 <= number < n"
  (let ((random-number (park-miller-randomizer)))
    ;; We subtract the randomly generated number from 1.0
    ;; before scaling so that we end up in the range
    ;; 0.0 <= x < 1.0, not 0.0 < x <= 1.0
    (* n (- 1.0d0 random-number))))

(defun random-integer (n)
  "Returns a pseudo-random integer in the range 0 ---> n-1."
  (let ((random-number (random-floating-point-number 1.0)))
    (floor (* n random-number))))

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

;;; Regression Problem for 0.5x**2

(defvar x)

(defun define-terminal-set-for-REGRESSION ()
 (values' (x :floating-point-random-constant))
)

(defun define-function-set-for-REGRESSION ()
  (values '(+ - * % )
          '(2 2 2 2 ))
)

(defun % (numerator denominator)
  "The Protected Division Function"
  (values (if (= 0 denominator) 1 (/ numerator denominator)))
)


(defun lg (base)
"The Protected Division Function"
 (values (if (= 0 base) 1 (log base)))
)

(defstruct REGRESSION-fitness-case
    independent-variable
    target
)

(defun define-fitness-cases-for-REGRESSION ()	;01
  (let (fitness-cases w x y this-fitness-case)	;02
  (setf w '(-5 -4 -3 -2 -1 0 1 2 3 4 5)) 
  (setf y '(25 16 9 4 1 0 1 4 9 16 25)) 
  (setf fitness-cases (make-array *number-of-fitness-cases*))	;03
    (format t "~%Fitness cases")	;04
    (dotimes (index *number-of-fitness-cases*)	;05
    (setf x (nth index  w))      
    (setf this-fitness-case (make-REGRESSION-fitness-case))	;07
      (setf (aref fitness-cases index) this-fitness-case)	;08
      (setf (REGRESSION-fitness-case-independent-variable	;09
             this-fitness-case)	;10
             x)	;11
      (setf (REGRESSION-fitness-case-target	;12
             this-fitness-case)	;13
            (nth index y))	;14
      (format t "~% ~D      ~D      ~D"	;15
            index	;16
            (float x)	;17
            (REGRESSION-fitness-case-target this-fitness-case))	;18
    )	;19
    (values fitness-cases)	;20
  )	;21
)	;22

(defun REGRESSION-wrapper (result-from-program)
  (values result-from-program)
)

(defun evaluate-standardized-fitness-for-REGRESSION	;01
           (program fitness-cases)	;02
  (let (raw-fitness hits standardized-fitness x target-value	;03
        difference value-from-program this-fitness-case)	;04
    (setf raw-fitness 0.0)	;05
    (setf hits 0)	;06
    (dotimes (index *number-of-fitness-cases*)	;07
      (setf this-fitness-case (aref fitness-cases index))	;08
      (setf   x	;09
            (REGRESSION-fitness-case-independent-variable	;10
                   this-fitness-case))	;11
      (setf target-value	;12
             (REGRESSION-fitness-case-target	;13
                   this-fitness-case))	;14
       (setf value-from-program	;15
             (REGRESSION-wrapper (eval program)))	;16
       (setf difference (abs (- target-value	;17
                                     value-from-program)))	;18
      (incf raw-fitness difference)	;19
      (when (< difference 0.01) (incf hits)))	;20
    (setf standardized-fitness raw-fitness)	;21
    (values standardized-fitness hits)	;22
  )	;23
)	;24

(defun define-parameters-for-REGRESSION ()
  (setf *number-of-fitness-cases* 11)
  (setf *max-depth-for-new-individuals* 2)
  (setf *max-depth-for-individuals-after-crossover* 1)
  (setf *fitness-proportionate-reproduction-fraction* 0.1)
  (setf *crossover-at-any-point-fraction* 0.2)
  (setf *crossover-at-function-point-fraction* 0.2)
  (setf *max-depth-for-new-subtrees-in-mutants* 1)
  (setf *method-of-selection* :tournament)
  (setf *method-of-generation* :ramped-half-and-half)
  (values)
)

(defun define-termination-criterion-for-REGRESSION	;01
       (current-generation	;02
        maximum-generations	;03
        best-standardized-fitness	;04
        best-hits)	;05
  (declare (ignore best-standardized-fitness))	;06
  (values	;07
    (or (>= current-generation maximum-generations)	;08
         (>= best-hits (* 2 *number-of-fitness-cases*)))	;09
  )	;10
)	;11

(defun REGRESSION ()
  (values 'define-function-set-for-REGRESSION
          'define-terminal-set-for-REGRESSION
          'define-fitness-cases-for-REGRESSION
          'evaluate-standardized-fitness-for-REGRESSION
          'define-parameters-for-REGRESSION
          'define-termination-criterion-for-REGRESSION
  )
)
(run-genetic-programming-system 'REGRESSION 1 400 2000)