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     (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 (declaim (ftype (function ((simple-array double-float (* *))
 207                            fixnum)
 208                           (values (array double-float (*)) &optional))
 209                 displace))
 210 (defun displace (simplex i)
 211   "Extract the i-th i=0..n point from the simplex."
 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 (declaim (ftype (function ((array double-float (*))
 219                            (array double-float (*))
 220                            double-float)
 221                           (values (array double-float (*)) &optional))
 222                 v-extrapolate))
 223 (defun v-extrapolate (vec-a vec-b alpha)
 224   "Extrapolate the two vectors like this: <1+alpha> a - alpha b."
 225   (v+ (v* vec-a (1+ alpha))
 226       (v* vec-b (- alpha))))
 227
 228 (declaim (ftype (function ((simple-array double-float (* *)) fixnum)
 229                           (values (simple-array double-float *) &optional))
 230                 calc-centroid))
 231 (defun calc-centroid (simplex worst)
 232   "Go through all n+1 points in the simplex except the worst and sum
 233 the coordinates. Return the centroid."
 234   (destructuring-bind (nr-points n)
 235       (array-dimensions simplex) ;; first dimension is bigger (n+1)
 236     (let ((result (make-array n :element-type 'double-float
 237                               :initial-element 0d0)))
 238       (dotimes (i nr-points)
 239         (unless (eq i worst)
 240           (v-inc result (displace simplex i))))
 241       (dotimes (j n)
 242         (setf (aref result j) (/ (aref result j) n)))
 243       (dformat t "~a~%" (list 'centroid result))
 244       result)))
 245
 246
 247 (declaim (ftype (function ((simple-array double-float (* *))
 248                            function
 249                            &key (:itmax fixnum)
 250                            (:ftol double-float)
 251                            (:temperature double-float))
 252                           (values double-float
 253                                   (simple-array double-float *)
 254                                   double-float
 255                                   double-float
 256                                   (simple-array double-float *)
 257                                   &optional))))
 258 (defun amoeba (simplex funk &key (itmax 500) (ftol 1d-5) (temperature 1000d0))
 259   (destructuring-bind (nr-points n)
 260       (array-dimensions simplex)
 261    (let* ((alpha 1d0) ;; reflected size
 262           (rho 2d0) ;;  expansion ratio
 263           (gamma .5d0) ;; contraction ratio
 264           (sigma .5d0) ;; shrink ratio
 265           (iteration 0)
 266           (best-ever (make-array n :element-type 'double-float))
 267           (best-ever-value 1d100)
 268           (values (make-array nr-points :element-type 'double-float)))
 269      ;; evaluate function on all vertices
 270      (dotimes (i nr-points)
 271        (setf (aref values i)
 272              (funcall funk (displace simplex i))))
 273      (tagbody
 274       label1
 275         (incf iteration)
 276         (multiple-value-bind (best best-value
 277                                    second-worst second-worst-value
 278                                    worst worst-value)
 279             (find-extrema values temperature)
 280           (declare (ignore second-worst))
 281           (when (< best-value best-ever-value) ;; store the best value
 282             (v-set best-ever (displace simplex best))
 283             (setf best-ever-value best-value))
 284           (dformat t "~a~%" (list 'best-ever best-ever best-ever-value))
 285           ;;      rtol=2 |h-l| / ( |h| + |l| + TINY)
 286           (let ((rtol (* 2d0 (/ (abs (- best-value worst-value))
 287                                 (+ (abs best-value)
 288                                    (abs worst-value)
 289                                    1d-20)))))
 290             (when (or (< rtol ftol)
 291                       (< itmax iteration))
 292               (return-from amoeba (values best-value
 293                                           (displace simplex best)
 294                                           rtol
 295                                           best-ever-value
 296                                           best-ever))))
 297
 298           (macrolet ((replace-worst (new-point new-value)
 299                        `(progn
 300                           (v-set (displace simplex worst) ,new-point)
 301                           (setf (aref values worst) ,new-value))))
 302             (let* ((centroid (calc-centroid simplex worst))
 303                    ;; reflect worst point around centroid
 304                    (reflected (v-extrapolate centroid (displace simplex worst) alpha))
 305                    (reflected-value (funcall funk reflected))
 306                    (reflected-heat (heat (- temperature) reflected-value)))
 307               (dformat t "~a~%" (list 'reflected reflected reflected-value))
 308               ;; if new point better than current best, expand the simplex
 309               (when (< reflected-heat best-value)
 310                 (let* ((expanded (v+ (v* reflected rho)
 311                                      (v* centroid (- 1 rho))))
 312                        (expanded-value (funcall funk expanded))
 313                        (expanded-heat (heat (- temperature) expanded-value)))
 314                   (dformat t "~a~%" (list 'expanded expanded expanded-value))
 315                   ;; keep whichever is best
 316                   (if (< expanded-heat reflected-value)
 317                       (replace-worst expanded expanded-value)
 318                       (replace-worst reflected reflected-value))
 319                   (go label1)))
 320               ;; new point lies between others, keep:
 321               (when (< reflected-heat second-worst-value)
 322                 (dformat t "keep~%")
 323                 (replace-worst reflected reflected-value)
 324                 (go label1))
 325               ;; new point is better than the worst one we had, contract:
 326               (when (< reflected-heat worst-value)
 327                 (let* ((contracted (v-extrapolate centroid
 328                                                   (displace simplex worst)
 329                                                   gamma))
 330                        (contracted-value (funcall funk contracted))
 331                        (contracted-heat (heat (- temperature) contracted-value)))
 332                   (dformat t "~a~%" (list 'contracted
 333                                          contracted contracted-value))
 334                   ;; use it if it helped
 335                   (when (< contracted-heat reflected-value)
 336                     (replace-worst contracted contracted-value)
 337                     (go label1))))
 338               ;; otherwise if reflected is worse than everything else or
 339               ;; contracted was worse than reflected shrink the whole
 340               ;; simplex around the best point
 341               (dotimes (i nr-points)
 342                 (unless (eq i best)
 343                   (v-set (displace simplex i)
 344                          (v+ (displace simplex best)
 345                              (v* (v- (displace simplex i)
 346                                      (displace simplex best))
 347                                  sigma)))
 348                   (setf (aref values i) (funcall funk (displace simplex i)))))
 349               (go label1))))))))
 350
 351 (declaim (ftype (function ((simple-array double-float (* *))
 352                            function
 353                            &key (:itmax fixnum) (:ftol double-float)
 354                            (:start-temperature double-float))
 355                           (values double-float
 356                                   (simple-array double-float *) &optional))))
 357 (defun anneal (simplex funk &key (itmax 500) (ftol 1d-5) (start-temperature 100d0))
 358   (let* ((m 30)
 359          (eps/m .02d0)
 360          (eps (* eps/m m))
 361          (temp start-temperature))
 362     (do ((count 0 (incf count)))
 363         (())
 364       (multiple-value-bind (value point rtol best-ever-value best-ever)
 365           (amoeba simplex funk :itmax m :ftol ftol
 366                   :temperature temp)
 367         (declare (ignore best-ever-value
 368                          best-ever))
 369         (when (or (< itmax count)
 370                   (< rtol ftol))
 371           (return-from anneal (values value point))))
 372       (setf temp (* (- 1 eps) temp)))))
 373
 374 ;; run the following code to test the downhill simplex optimizer on a
 375 ;; 2d function:
 376 #+nil
 377 (time (let ((start (make-array 2 :element-type 'double-float
 378                           :initial-contents (list 1.5d0 1.5d0))))
 379    (anneal (make-simplex start 1d0)
 380            #'rosenbrock
 381            :ftol 1d-5)))
 382
 383 (declaim (ftype (function ((simple-array double-float *)
 384                            double-float)
 385                           (values (simple-array double-float (* *)) &optional))
 386                 make-simplex))
 387 (defun make-simplex (start step)
 388   (let* ((n (array-dimension start 0))
 389          (result (make-array (list (1+ n) n) :element-type 'double-float)))
 390     (dotimes (j (1+ n))
 391       (v-set (displace result j) start))
 392     ;; change the first n points by step
 393     ;; the last point in result stays unaltered
 394     (dotimes (j n)
 395       (incf (aref result j j) step))
 396
 397     result))
 398
 399 ;; (declaim (ftype (function (double-float)
 400 ;;                        (values double-float &optional))
 401 ;;              sq))
 402 ;; (defun sq (x)
 403 ;;   (* x x))
 404
 405 ;; (declaim (ftype (function ((array double-float *))
 406 ;;                        (values double-float &optional))
 407 ;;              rosenbrock))
 408 ;; (defun rosenbrock (p)
 409 ;;   (let* ((x (aref p 0))
 410 ;;       (y (aref p 1))
 411 ;;       (result (+ (sq (- 1 x))
 412 ;;                  (* 100 (sq (- y (sq x)))))))
 413 ;;     (dformat t "~a~%" (list 'rosenbrock p result))
 414 ;;     result))
 415 ;; #+nil
 416 ;; (rosenbrock (make-array 2 :element-type 'double-float
 417 ;;                       :initial-contents (list 1.5d0 1.5d0)))