simplex-anneal.lisp

   1 (defpackage :simplex-anneal
   2   (:use :cl)
   3   (:export #:make-simplex
   4            #:anneal))
   5 (in-package :simplex-anneal)
   6
   7 ;; nr3.pdf p.526 10.5 downhill simplex
   8
   9 #|| A simplex in N dimensions has N+1 points, i.e. a triangle in 2D.
  10 You give one ND point. The algorithm constructs the simplex around
  11 it (distances depend on the problems characteristic length scale)
  12  and goes downhill.
  13
  14 There are five different types of steps:
  15 a) reflection (conserving the volume)
  16 b) reflection and expansion
  17 c) contraction
  18 d) multiple-contraction
  19
  20 We keep track of highest and next to highest
  21
  22 Termination:
  23 Step is fractionally smaller than some tolerance tol.
  24 OR
  25 Decrease in function value is fractionally smaller than some tolerance tol.
  26
  27 tol > sqrt(eps), tol can be \approx eps
  28
  29 If a minimum is found, restart the routine with 2 new points and the
  30 current best one.
  31
  32 3x2 matrix defines simplex in 2D (rows are coordinates).
  33
  34 Store 3 function values.
  35 Determine highest, next highest and lowest (best) point.
  36
  37 Fractional range from lowest value l to highest value h:
  38 rtol=2 |h-l| / ( |h| + |l| + TINY)
  39 TINY is 1e-20 in cvd
  40
  41 Extrapolate by factor -1 through the face oposite of the high point.
  42  - If result better than best point extrapolate by two in same direction.
  43  - Else if worse than second highest do 1d contraction of 0.5,
  44    searching for lower one
  45
  46 amotry does the extrapolation
  47 it gets the points psum, the values, a double fac
  48
  49 fac1=(1-fac)/n
  50 fac2=fac1-fac   = 1/n - fac/n -fac = 1/n - fac ( 1/n - 1 )
  51
  52 mirrored point
  53 j goes over n:
  54 try_j= psum_j *fac1 -p [high]_j *fac2
  55
  56 if the new value is better than the worst (highest) replace that point
  57
  58 get_psum has 2 values for 2D and contains the sum of the coordinates
  59 of all points
  60
  61 in cvd the reflection is around the centroid
  62 alpha 1  reflected size
  63 rho 2  expansion ratio
  64 gamma .5 contraction ratio
  65 sigma .5 shrink ratio
  66 epsilon sqrt(eps) tolerance
  67
  68 For simulated annealing add a positive logarithmic distributed random
  69 variable proportional to the temperature: (* temp (log (random))) to
  70 every stored function value of the vertices. A similar variable is
  71 subtracted from every point that is tried.
  72
  73 At finite T the simplex has approximately the shape of the region,
  74 that can be reached at this temperature. It tumbles with Brownian
  75 motion within that region.
  76
  77 Possible anneling temperature schedules:
  78 - T=(1-eps)*T after every m moves. eps/m is determined by experiment.
  79 - Budget K moves. Reduce T after m moves to T_0*(1-k/K)^alpha with k the
  80   count of moves so far and a constant alpha=1,2 or 4 or something. alpha
  81   depends on the statistical distribution of relative minima at
  82   various depths. Larger alpha spend more iterations at low temperature.
  83 - After every m moves set T to beta*(f_1-f_b) with beta\approx 1, f_1 is the
  84   smallest function value currently in the simplex and f_b is the best
  85   value ever encountered. But never reduce T by more than some
  86   fraction gamma.
  87
  88 An occasional restart may help.
  89
  90 ||#
  91
  92 (declaim (optimize (speed 2) (safety 3) (debug 3)))
  93
  94 (defparameter *print-debug-message* nil)
  95
  96 (defmacro dformat (&rest rest)
  97   `(progn ;; when simplex-anneal::*print-debug-message*
  98      (format ,@rest)))
  99
 100 #+nil
 101 (dformat t "~a~%" (list 'test))
 102
 103 (declaim (ftype (function (double-float double-float)
 104                           (values double-float &optional))
 105                 heat))
 106 (defun heat (temperature value)
 107   (- value (* temperature (log (random 1d0)))))
 108
 109 (defmacro aref-heat (temperature array index)
 110   "Access an array and change values with random fluctation."
 111   `(heat ,temperature (aref ,array ,index)))
 112
 113 (declaim (ftype (function ((simple-array double-float *) double-float)
 114                           (values fixnum double-float
 115                                   fixnum double-float
 116                                   fixnum double-float &optional))
 117                 find-extrema))
 118 (defun find-extrema (values temperature)
 119   "Find the best, second-worst and worst point in the array of simplex
 120 function values."
 121   (let* ((best 0)
 122          (worst best)
 123          (bestv (aref-heat temperature values best))
 124          (worstv bestv)
 125          (second-worstv 0d0)
 126          (second-worst 0))
 127     ;; be careful to maintain the right relationship between worst
 128     ;; and second-worst when initializing: if 0 could be the worst
 129     ;; one use something better, so check 1 and use this if its
 130     ;; better.
 131     (let ((help (aref-heat temperature values 1)))
 132       (if (< help worstv)
 133           (setf second-worst 1 ;; 1 is better
 134                 second-worstv help)
 135           (setf second-worst 0 ;; 1 is worse
 136                 second-worstv worstv
 137                 worst 1
 138                 worstv help)))
 139     ;; we don't have to loop over 0, so start with 1
 140     (loop for i from 1 below (length values) do
 141       (let ((v (aref-heat temperature values i)))
 142         (when (< v bestv)
 143           (setf bestv v
 144                 best i))
 145         (when (< worstv v)
 146           (setf worstv v
 147                 worst i))
 148         (unless (eq i worst)
 149           (when (< second-worstv v)
 150             (setf second-worstv v
 151                   second-worst i)))))
 152 #+nil    (dformat t "~a~%" (list
 153                        'best best bestv
 154                        'second-worst second-worst second-worstv
 155                        'worst worst worstv))
 156     (values best bestv second-worst second-worstv worst worstv)))
 157
 158 (declaim (ftype (function ((array double-float *)
 159                            (array double-float *))
 160                           (values (simple-array double-float *) &optional))
 161                 v+ v- v-inc))
 162 (defun v+ (a b)
 163   (let* ((n (array-dimension a 0))
 164          (result (make-array n :element-type 'double-float)))
 165     (dotimes (i n)
 166       (setf (aref result i) (+ (aref a i) (aref b i))))
 167     result))
 168
 169 (defun v- (a b)
 170   (let* ((n (array-dimension a 0))
 171          (result (make-array n :element-type 'double-float)))
 172     (dotimes (i n)
 173       (setf (aref result i) (- (aref a i) (aref b i))))
 174     result))
 175
 176 (defun v-inc (vec-result vec)
 177   "Increase the values in VEC-RESULT with values from VEC. Return the
 178 modified VEC-RESULT."
 179   (let* ((n (array-dimension vec-result 0)))
 180     (dotimes (i n)
 181       (incf (aref vec-result i) (aref vec i)))
 182     vec-result))
 183
 184 (declaim (ftype (function ((array double-float *)
 185                            (array double-float *))
 186                           (values (array double-float *) &optional))
 187                 v-set))
 188 (defun v-set (vec-result vec)
 189   "Set the values in the vector VEC-RESULT to the contents of VEC."
 190   (let* ((n (array-dimension vec-result 0)))
 191     (dotimes (i n)
 192       (setf (aref vec-result i) (aref vec i)))
 193     vec-result))
 194
 195 (declaim (ftype (function ((array double-float *)
 196                            double-float)
 197                           (values (simple-array double-float *) &optional))
 198                 v*))
 199 (defun v* (vec scalar)
 200   (let* ((n (array-dimension vec 0))
 201          (result (make-array n :element-type 'double-float)))
 202     (dotimes (i n)
 203       (setf (aref result i) (* scalar (aref vec i))))
 204     result))
 205
 206 (defun displace-reference (simplex i)
 207   "Extract the i-th i=0..n point from the simplex. Changing the result
 208 will modify the simplex as well."
 209   (declare ((simple-array double-float 2) simplex)
 210            (fixnum i)
 211            (values (array double-float 1) &optional))
 212   (let ((n (array-dimension simplex 1))) ;; we want the small index
 213     (make-array n
 214                 :element-type 'double-float
 215                 :displaced-to simplex
 216                 :displaced-index-offset (* i n))))
 217
 218 (defun displace (simplex i)
 219   "Extract the i-th i=0..n point from the simplex."
 220   (declare ((simple-array double-float 2) simplex)
 221            (fixnum i)
 222            (values (simple-array double-float 1) &optional))
 223   (let ((n (array-dimension simplex 1))) ;; we want the small index
 224     ;; force copy into simple-array
 225     (make-array n :element-type 'double-float
 226                 :initial-contents
 227                 (displace-reference simplex i))))
 228
 229
 230
 231 (declaim (ftype (function ((array double-float (*))
 232                            (array double-float (*))
 233                            double-float)
 234                           (values (array double-float (*)) &optional))
 235                 v-extrapolate))
 236 (defun v-extrapolate (vec-a vec-b alpha)
 237   "Extrapolate the two vectors like this: <1+alpha> a - alpha b."
 238   (v+ (v* vec-a (1+ alpha))
 239       (v* vec-b (- alpha))))
 240
 241 (declaim (ftype (function ((simple-array double-float (* *)) fixnum)
 242                           (values (simple-array double-float *) &optional))
 243                 calc-centroid))
 244 (defun calc-centroid (simplex worst)
 245   "Go through all n+1 points in the simplex except the worst and sum
 246 the coordinates. Return the centroid."
 247   (destructuring-bind (nr-points n)
 248       (array-dimensions simplex) ;; first dimension is bigger (n+1)
 249     (let ((result (make-array n :element-type 'double-float
 250                               :initial-element 0d0)))
 251       (dotimes (i nr-points)
 252         (unless (eq i worst)
 253           (v-inc result (displace simplex i))))
 254       (dotimes (j n)
 255         (setf (aref result j) (/ (aref result j) n)))
 256 ;;      (dformat t "~a~%" (list 'centroid result))
 257       result)))
 258
 259 ;; 1/n! * det(v1-v0,v2-v0,...,vn-v0)
 260 ;; only works for 2d coordinates A=((a b)(c d)), det A = ad-bc:
 261 ;; a = v1_x - v0_x
 262 (defun simplex-volume (simplex)
 263   (declare ((simple-array double-float (3 2)) simplex)
 264            (values double-float &optional))
 265   (let ((a (- (aref simplex 1 0) (aref simplex 0 0)))
 266         (c (- (aref simplex 1 1) (aref simplex 0 1)))
 267         (b (- (aref simplex 2 0) (aref simplex 0 0)))
 268         (d (- (aref simplex 2 1) (aref simplex 0 1))))
 269     (- (* a d) (* b c))))
 270
 271
 272 (defun amoeba (simplex funk &key (itmax 500) (ftol 1d-5) (temperature 1000d0)
 273                (simplex-min-size .1d0))
 274   (declare ((simple-array double-float 2) simplex)
 275            (function funk)
 276            (fixnum itmax)
 277            (double-float ftol temperature simplex-min-size)
 278            (values double-float
 279                    (simple-array double-float 1)
 280                    double-float
 281                    double-float
 282                    (simple-array double-float 1)
 283                    &optional))
 284   (destructuring-bind (nr-points n)
 285       (array-dimensions simplex)
 286    (let* ((alpha 1d0) ;; reflected size
 287           (rho 2d0) ;;  expansion ratio
 288           (gamma .5d0) ;; contraction ratio
 289           (sigma .5d0) ;; shrink ratio
 290           (iteration 0)
 291         ;;  (min-volume (* simplex-min-size simplex-min-size simplex-min-size))
 292           (best-ever (make-array n :element-type 'double-float))
 293           (best-ever-value 1d100)
 294           (vals (make-array nr-points :element-type 'double-float)))
 295      ;; evaluate function on all vertices
 296      (dotimes (i nr-points)
 297        (setf (aref vals i)
 298              (funcall funk (displace simplex i))))
 299      (tagbody
 300       label1
 301         (incf iteration)
 302         (multiple-value-bind (best best-value
 303                                    second-worst second-worst-value
 304                                    worst worst-value)
 305             (find-extrema vals temperature)
 306           (declare (ignore second-worst))
 307           ;;(format t "~f~%" best-value)
 308           (when (< best-value best-ever-value) ;; store the best value
 309             (v-set best-ever (displace-reference simplex best))
 310             (setf best-ever-value best-value))
 311 ;;        (dformat t "~a~%" (list 'best-ever best-ever best-ever-value))
 312           ;;      rtol=2 |h-l| / ( |h| + |l| + TINY)
 313           ;; fractional range from highest to lowest and return if satisfied
 314           (let ((rtol (* 2d0 (/ (abs (- best-value worst-value))
 315                                 (+ (abs best-value)
 316                                    (abs worst-value)
 317                                    1d-20)))))
 318             (when (or (< rtol ftol)
 319                       ;; should only be checked for low temperatures
 320                       ;;(< (abs (simplex-volume simplex)) min-volume)
 321                       (< itmax iteration))
 322               (return-from amoeba (values best-value
 323                                           (displace simplex best)
 324                                           rtol
 325                                           best-ever-value
 326                                           best-ever))))
 327
 328           (macrolet ((replace-worst (new-point new-value)
 329                        `(progn
 330                           (v-set (displace-reference simplex worst) ,new-point)
 331                           (setf (aref vals worst) ,new-value))))
 332             (let* ((centroid (calc-centroid simplex worst))
 333                    ;; reflect worst point around centroid
 334                    (reflected (v-extrapolate centroid (displace simplex worst) alpha))
 335                    (reflected-value (funcall funk reflected))
 336                    (reflected-heat (heat (- temperature) reflected-value)))
 337         ;;      (dformat t "~a~%" (list 'reflected reflected reflected-value))
 338               ;; if new point better than current best, expand the simplex
 339               (when (< reflected-heat best-value)
 340                 (let* ((expanded (v+ (v* reflected rho)
 341                                      (v* centroid (- 1 rho))))
 342                        (expanded-value (funcall funk expanded))
 343                        (expanded-heat (heat (- temperature) expanded-value)))
 344                 ;;  (dformat t "~a~%" (list 'expanded expanded expanded-value))
 345                   ;; keep whichever is best
 346                   (if (< expanded-heat reflected-value)
 347                       (replace-worst expanded expanded-value)
 348                       (replace-worst reflected reflected-value))
 349                   (go label1)))
 350               ;; new point lies between others, keep:
 351               (when (< reflected-heat second-worst-value)
 352         ;;      (dformat t "keep~%")
 353                 (replace-worst reflected reflected-value)
 354                 (go label1))
 355               ;; new point is better than the worst one we had, contract:
 356               (when (< reflected-heat worst-value)
 357                 (let* ((contracted (v-extrapolate centroid
 358                                                   (displace simplex worst)
 359                                                   gamma))
 360                        (contracted-value (funcall funk contracted))
 361                        (contracted-heat (heat (- temperature) contracted-value)))
 362                 #+nil  (dformat t "~a~%" (list 'contracted
 363                                          contracted contracted-value))
 364                   ;; use it if it helped
 365                   (when (< contracted-heat reflected-value)
 366                     (replace-worst contracted contracted-value)
 367                     (go label1))))
 368               ;; otherwise if reflected is worse than everything else or
 369               ;; contracted was worse than reflected shrink the whole
 370               ;; simplex around the best point
 371               (dotimes (i nr-points)
 372                 (unless (eq i best)
 373                   (v-set (displace-reference simplex i)
 374                          (v+ (displace simplex best)
 375                              (v* (v- (displace simplex i)
 376                                      (displace simplex best))
 377                                  sigma)))
 378                   (setf (aref vals i) (funcall funk (displace simplex i)))))
 379               (go label1))))))))
 380
 381 (defun anneal (simplex funk &key (itmax 500) (ftol 1d-5) (start-temperature 100d0)
 382                (eps/m .02d0) (simplex-min-size .1d0))
 383   (declare ((simple-array double-float 2) simplex)
 384            (function funk)
 385            (fixnum itmax)
 386            (double-float ftol start-temperature eps/m simplex-min-size)
 387            (values double-float (simple-array double-float 1) &optional))
 388   (let* ((m 30)
 389          (eps (* eps/m m))
 390          (temp start-temperature))
 391     (do ((count 0 (incf count)))
 392         (())
 393       (multiple-value-bind (value point rtol best-ever-value best-ever)
 394           (amoeba simplex funk :itmax m :ftol ftol
 395                   :temperature temp :simplex-min-size simplex-min-size)
 396         (declare (ignore best-ever-value
 397                          best-ever))
 398         (when (or (< itmax count)
 399                   (< rtol ftol))
 400           (return-from anneal (values value point))))
 401       (setf temp (* (- 1 eps) temp)))))
 402
 403 ;; run the following code to test the downhill simplex optimizer on a
 404 ;; 2d function:
 405 #+nil
 406 (time
 407  (progn
 408    (defun sq (x)
 409      (declare (double-float x)
 410               (values double-float &optional))
 411      (* x x))
 412    (defun rosenbrock (p)
 413      (declare ((array double-float 1) p)
 414               (values double-float &optional))
 415      (let* ((x (aref p 0))
 416             (y (aref p 1))
 417             (result (+ (sq (- 1 x))
 418                        (* 100 (sq (- y (sq x)))))))
 419        #+nil (format t "~a~%" (list 'rosenbrock p result))
 420      result))
 421    (with-open-file (*standard-output* "/dev/shm/a"
 422                                       :if-exists :supersede
 423                                       :if-does-not-exist :create
 424                                       :direction :output)
 425     (let ((start (make-array 2 :element-type 'double-float
 426                              :initial-contents (list 1.5d0 1.5d0))))
 427       (anneal (make-simplex start 1d0)
 428               #'rosenbrock
 429               :start-temperature 100d0
 430               :itmax 100
 431               :ftol 1d-3
 432               :eps/m .01d0)))))
 433
 434 #+nil
 435 (rosenbrock (make-array 2 :element-type 'double-float
 436                         :initial-contents (list 1.5d0 1.5d0)))
 437
 438
 439 (defun make-simplex (start step)
 440   (declare ((simple-array double-float 1) start)
 441            (double-float step)
 442            (values (simple-array double-float 2) &optional))
 443   (let* ((n (array-dimension start 0))
 444          (result (make-array (list (1+ n) n) :element-type 'double-float)))
 445     (dotimes (j (1+ n))
 446       (v-set (displace-reference result j) start))
 447     ;; change the first n points by step
 448     ;; the last point in result stays unaltered
 449     (dotimes (j n)
 450       (incf (aref result j j) step))
 451     result))