Differential evolution in LISP

de.lisp

#| 
############################################
de.lisp: Differential evolution
Tim Menzies, [email protected], http://menzies.us
(C) 2012 GPL 3.0

To use: (load "de.lisp")

For more on differential evolution, 
see http://tinyurl.com/83d9daa
############################################### 
|#

(defun de (&key 
	   (cf 0.5) (f 0.7) (np 10000)  
	   (d (make-d :cf cf :f f :np np)))
    "CL-USER> (de :f 0.25)
    100000 [100%] 22%  0.52 (  6  5  3  7  3  4  4  4  5  3)
     99998 [100%] 38%  0.72 (  9  7  5  9  7  5  4  6  9  7)
     99986 [100%]  4%  0.75 (  4  8  8  5  8  8  9  8  5 10)
     99520 [100%]  0%  0.75 (  5  9  7  8  5  7 10  5  7  9)
     99502 [100%]  1%  0.76 (  7 10  7 10  6  8  8  7  6  4)
     99481 [ 99%]  0%  0.76 (  9 10  6  6  6  5  9 10  7  6)
     99193 [ 99%]  3%  0.78 (  9  8 10  7  9  4  7  9  7  6)
     98706 [ 99%]  3%  0.80 (  7  6 10  9  6 10  8  6 10  7)
     97091 [ 97%]  2%  0.82 (  9  7  6  9  9  6  8  8 10  6)
     96009 [ 96%]  4%  0.85 (  9  6  9  8  7 10  8  8  7  8)
     80994 [ 81%]  3%  0.88 ( 10  7 10  8  8  7 10 10 10  9)
     66706 [ 67%]  2%  0.90 (  9  7  9 10 10 10  8  8  8 10)
     50099 [ 50%]  1%  0.90 (  9  8 10  8  8 10  7  8  9  9)
     48810 [ 49%]  2%  0.92 ( 10 10  8  9 10  9  8  8  9 10)
     20360 [ 20%]  2%  0.94 ( 10  9  9  9  9  8 10  9  9  8)
    want : (9 9 9 9 9 9 9 9 9 9)
    got  : (10 9 9 9 9 8 10 9 9 8)
    NIL
    CL-USER> 
    ; No value
    CL-USER> (de :f 0.5)
    100000 [100%] 22%  0.52 (  5  5  3  7  3  4  4  4  6  3)
     99998 [100%] 36%  0.71 (  9  7  4  9  8  6  4  6  8  6)
     99986 [100%]  4%  0.74 (  4  7 10  6  8  8  9 10  4 10)
     99913 [100%]  1%  0.75 (  6  6  6  5  9  5  9  9  7  9)
     99520 [100%]  1%  0.75 (  6 10  7  8  5  7 10  5  6  9)
     99481 [ 99%]  5%  0.79 (  8  9  6  7  6  6  9 10  8  6)
     98706 [ 99%]  1%  0.80 (  7  6 10  9  6 10  8  6 10  7)
     97091 [ 97%]  2%  0.82 (  9  7  6  9  9  6  8  8 10  6)
     96009 [ 96%]  2%  0.83 (  9  6  9  8  7 10  8  6  7  9)
     91622 [ 92%]  2%  0.85 (  8  9  8 10  8  6  6  9 10  8)
     87396 [ 87%]  1%  0.86 (  5  9 10  9  9  8 10  9  8  8)
     86596 [ 87%]  3%  0.88 (  9  9  8  9 10  6 10 10  8 10)
     71702 [ 72%]  1%  0.89 (  7 10 10  8 10 10  8 10 10 10)
     63911 [ 64%]  1%  0.89 ( 10 10  7  9 10 10 10  8  8 10)
     62597 [ 63%]  2%  0.91 ( 10  8  9  8 10 10  8  8  8  8)
     36319 [ 36%]  1%  0.92 ( 10 10  9  8 10 10  9 10 10  9)
     19782 [ 20%]  1%  0.92 (  9  9  8  8 10  8  8  9 10  9)
    want : (9 9 9 9 9 9 9 9 9 9)
    got  : (9 9 8 8 10 8 8 9 10 9)
    NIL
    CL-USER> (de :f 0.75)
    100000 [100%] 20%  0.51 (  4  5  3  7  3  4  4  4  6  3)
     99998 [100%] 30%  0.67 (  9  7  3  9 10  6  4  6  8  4)
     99991 [100%]  5%  0.70 ( 10  8  7  4  6  4  6  6  8  7)
     99986 [100%]  5%  0.74 (  4  6 10  7  8  8  9 10  4 10)
     99752 [100%]  1%  0.75 (  9  9  5  9  4  7  7  8  6  7)
     99520 [100%]  3%  0.77 (  8 10  7  8  5  7 10  5  6  9)
     99481 [ 99%]  4%  0.80 (  8  8  6  8  6  6  9 10  9  6)
     98706 [ 99%]  0%  0.80 (  7  6 10  9  6 10  8  6 10  7)
     97091 [ 97%]  2%  0.82 (  9  7  6  9  9  6  8  8 10  6)
     91622 [ 92%]  2%  0.84 (  8  9  7 10  8  6  6  9 10  8)
     91282 [ 91%]  0%  0.84 (  7  6  7 10 10  8  8  9  7 10)
     90030 [ 90%]  0%  0.84 ( 10  6  8 10 10 10 10  6  9 10)
     89737 [ 90%]  3%  0.87 (  7 10  7  8  7  8 10 10  9  8)
     86596 [ 87%]  1%  0.88 (  8 10  7  9 10  7 10 10  8 10)
     77255 [ 77%]  1%  0.89 ( 10  7  9  7 10 10  8  9  9  9)
     74050 [ 74%]  1%  0.90 (  9 10 10  8 10  9 10  8  8  7)
     67554 [ 68%]  1%  0.91 ( 10 10  9  7 10  9  8  8  9  9)
     48810 [ 49%]  1%  0.92 (  9 10  9 10 10  9  8 10  8 10)
     17140 [ 17%]  1%  0.92 ( 10 10  8  9 10 10  9  9  9 10)
    want : (9 9 9 9 9 9 9 9 9 9)
    got  : (10 10 8 9 10 10 9 9 9 10)
    "
  (reset-seed)
  (dotimes (i (d-np d)) 
    (new d))
  (let* ((zero (any d))
	 (runs (1+ (* 10 (d-np d))))
	 (got  (de-run runs runs
		       zero 
		       (closest zero d) 
		       d)))
    (format t "~&want : ~a~%got  : ~a~%"  
	    (d-goal d) got)))

(defun de-run (n runs best best-score d)
  "run n times"
  (if (zerop n)
      best
      (de-step (- n 1) runs best best-score d)))

(defun de-step (n runs best best-score d)
  "for each member of the population, perhaps
    replace it with a mutated child"
  (let* ((id           (mod n (d-np d)))
	 (parent       (gethash id (d-all d) (any d)))
	 (child        ;(any d) 
		       (candidate d parent)
	               )
	 (parent-score (closest parent d))
	 (child-score  (closest child  d))
	 (winner       parent)
	 (winner-score parent-score))
    (when (> child-score parent-score)
      (setf winner                 child
	    winner-score           child-score
	    (gethash id (d-all d)) child))
    (if (<= winner-score best-score)
	(de-run n runs best   best-score   d)
	(progn
	  (format t "~&~6d [~3d%]~3d% ~5,2f (~{~3d~})~%" 
		  n
		  (round (* 100 (/ n runs)))
		  (round 
		   (* 100 (/ (- winner-score best-score)
			     best-score)))
		  winner-score winner)
	  (de-run n runs winner winner-score d)))))

(defun candidate (d parent)
  "The DE mutator candidate = any + f*(another - yetAnother)"
  (labels ((candidate1 (x y z r &aux (f (d-f d)))
	     (trim r (round (+ x (* f (- y z))))))
	   (cross-over (parent c)
	     (if (<= (randf 1.0) (d-cf d)) parent c)))
    (mapcar #'cross-over parent 
	    (mapcar #'candidate1 
		    (any d) (any d) (any d) (d-dd d)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defstruct range 
 " defining ranges for each attribute"
  (name (gensym "NAME"))
  (w 1)
  (min 0) 
  (max 10))

(defmethod any ((r range))
  "pull a random number from this range"
  (with-slots (min max) r
    (round (+ min (* (- max min) (randf 1.0))))))

(defmethod trim ((r range) n)
  "demand that a number falls within range"
  (with-slots (min max) r
    (cond ((> n max) max)
	  ((< n min) min)
	  (t         n))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defstruct d 
  "Defining the working memory for our de"
  (cf 0.5)
  (f  0.3)
  (np 10000)
  (dd (loop for x from 1 to 10
	 collect (make-range  :min 0 :max 10)))
  (goal '(9 9 9 9 9 9 9 9 9 9)) ; <=== secret goal
  (all   (make-hash-table)))


(defmethod any-id ((d d))
  "Select any id of the population"
  (1+ (randi  (hash-table-count (d-all d)))))

(defmethod  any ((d d))
  "Select any member of the population"
  (or (gethash (any-id d) (d-all d))
      (any d)))

(defmethod new ((x d))
  "Create a new member of the population,store 
   at old size + 1"
  (with-slots (dd n all) x
    (let ((where (+ 1 (hash-table-count all))))
      (setf (gethash where all)
	    (mapcar #'any dd)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun closest (x d)
  "objective function to be maximizes:
    compliment of distance between x and goal"
  (- 1 (dist x (d-goal d) (d-dd d))))

(defun !closest (&aux (d (make-d)))
  "testing closest. should return 1"
  (closest (d-goal d) d))

(defun dist(x y dd &aux (sum 0) (ws 0))
  "normalized distance between two things"
  (labels ((norm (n r)
	     (/ (- n (range-min r)) 
		(- (range-max r) (range-min r))))
	   (inc (a b r)
	     (let ((w (range-w r)))
	       (incf ws w)
	       (incf sum 
		     (* w (expt (- (norm a r) 
				   (norm b r)) 
				2))))))
    (mapc #'inc x y dd)
    (/ (sqrt sum) 
       (sqrt ws))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; utils. a random number generator where you
; can control the seed

(let* ((seed0      10013)
       (seed       seed0)
       (multiplier 16807.0d0)
       (modulus    2147483647.0d0))
  (defun reset-seed ()  (setf seed seed0))
  (defun randf      (n) (* n (- 1.0d0 (park-miller-randomizer))))
  (defun randi      (n) (floor (* n (/ (randf 1000.0) 1000))))
  (defun park-miller-randomizer ()
    "cycle= 2,147,483,646 numbers"
    (setf seed (mod (* multiplier seed) modulus))
    (/ seed modulus))
)

    

To see the effects of turning off different parts of this algorithm, in DE-STEP:

• I replaced (CANDIDATE D PARENT)
• with (ANY D).

This means that we no longer explore new solutions that are some interpolation of old solutions (instead, we just try combining all the existing stuff).

In this variant, we never rise above 83% of max score.

Conclusion: interpolation is useful.

CL-USER> (de)
100000 [100%]  8%  0.46 (  2  5  3  7  3  4  4  4  3  3)
 99998 [100%] 22%  0.57 (  9  7  1  9 10  3  4  6 10  2)
 99994 [100%] 22%  0.69 (  5  8  9  9  5  9  9  2  7  6)
 99991 [100%]  2%  0.70 ( 10  8  7  4  6  4  6  6  8  7)
 99946 [100%]  1%  0.71 (  8  8  7  8 10  5  3  7  6  6)
 99752 [100%]  5%  0.75 (  9  9  5  9  4  7  7  8  6  7)
 99340 [ 99%]  6%  0.80 (  9  8  8  8  7  9  9 10  6  4)
 99337 [ 99%]  4%  0.83 (  6 10  9  7  7  8  9  6  8  8)
want : (9 9 9 9 9 9 9 9 9 9)
got  : (6 10 9 7 7 8 9 6 8 8)

In another variant, I stopped replacing the parents with the candidate. In this variant, we also did not get very close to the max score.

Conclusion: pruning is useful.

CL-USER> (de)
100000 [100%] 20%  0.51 (  4  5  3  7  3  4  4  4  6  3)
 99998 [100%] 27%  0.65 (  9  7  3  9  6  3  4  6 10  8)
 99991 [100%]  8%  0.70 ( 10  8  7  4  6  4  6  6  8  7)
 99946 [100%]  1%  0.71 (  8  8  7  8 10  5  3  7  6  6)
 99819 [100%]  6%  0.76 (  6 10  8  5  7  5  6  9  7  9)
 99568 [100%]  1%  0.76 (  8  9  4  6  7  5  8  9 10  9)
 99193 [ 99%]  3%  0.78 (  9  8 10  7  9  4  7  9  7  6)
 98706 [ 99%]  3%  0.80 (  7  6 10  9  6 10  8  6 10  7)
 97928 [ 98%]  6%  0.85 (  8  8 10  7  7 10  7  7  8  8)
 86452 [ 86%]  0%  0.86 (  6  9  7  8  9  7 10  9 10  8)
 85436 [ 85%]  1%  0.87 (  9  8  7 10  8  9  9  8  6  8)
 30802 [ 31%]  1%  0.87 (  9 10  8  8 10  8 10  6 10  9)
want : (9 9 9 9 9 9 9 9 9 9)
got  : (9 10 8 8 10 8 10 6 10 9)

Finally, I played with the size of the population. The above examples used :np 10000. In the following, I tried 1% of that, i.e. :np 100. Note that we achieve nearly all what can be achieved above, in 1% of the time.

CL-USER> (de :np 100)
  1000 [100%] 40%  0.61 ( 10  1  5  8  7  1  8  9  9  7)
   997 [100%]  4%  0.63 ( 10  6  3  4  6  5  9 10  3  7)
   993 [ 99%]  4%  0.66 (  4  3  9  4  6 10  5  9  8  7)
   988 [ 99%]  5%  0.69 (  8  6  5  6 10  5  6  8  6  4)
   968 [ 97%]  5%  0.73 (  9  6 10 10  6  3  8  9  8  5)
   959 [ 96%]  1%  0.74 (  8  5  6  3  7  9 10 10 10 10)
   913 [ 91%]  5%  0.77 (  9  5 10  7  9  6  8  8  7  5)
   871 [ 87%]  7%  0.83 ( 10  9 10  6  9  6 10  7  7  9)
   473 [ 47%]  4%  0.87 (  9 10  9  6 10  7 10 10  9 10)
   407 [ 41%]  0%  0.87 (  7 10  9  6 10 10  9 10  9  9)
   234 [ 23%]  1%  0.88 ( 10  9 10  9 10  6  9 10 10 10)
   121 [ 12%]  3%  0.90 ( 10 10 10 10 10 10  9  9  7  9)
    30 [  3%]  2%  0.92 (  9  9  9  9  9 10 10 10  9  7)
want : (9 9 9 9 9 9 9 9 9 9)
got  : (9 9 9 9 9 10 10 10 9 7)

Should really move all and n and best and best-score into a runtime wme that is carried around de-run (fewer arguments need passing abound).

Should add more control variables so, on the command-line can:

• turn interpolation off/on
• turn replacing the parent off/on

Should really declare the static a global *d*. ease of reference

instances should be altered to be a structure that separates decision space from objective space e.g.

#S(eg :dec (1 2 3 4 5) :obj (23 42))

and the obj values are set when the instance is created. (so we never have to explicitly score the example with de-step)

Also, certain domain specific stuff needs to be isolated:

• definition of the dd
• objective function

Needs to be extended to handle multi-objective reasoning (as per Robic and Filipic. Which means implementing their three variants:

• DEMO/parent : already done. See my code.
• DEMO/closest/dec. Works in the same way as DEMO/parent, with the exception that the candidate solution is compared to the most similar individual in decision space. If it dominates it, the candidate replaces this individual, otherwise it is treated in the same way as in DEMO/parent. The applied similarity measure is the Euclidean distance between the two solutions in decision space.
• DEMO/closet/obj: In the third variant, DEMO/closest/obj, the candidate is compared to the most similar individual in objective space.

Note the runtime hit for the last two (searching through all others). We could reduce that, possibly, by augmented the most similar individual function to be most similar individual in a random sample of X members of decision/objective space.

Finally, the standard de cross-over function demands that at least one attribute from the parent (chosen at random) exists in the candidate.