From df8cdd72f63e269140e193017ba52b644db5fbe0 Mon Sep 17 00:00:00 2001 From: Robert Dodier Date: Thu, 28 Jul 2022 09:11:36 -0700 Subject: [PATCH] Commit new share package nelder_mead: Nelder-Mead algorithm for minimization without derivatives. Common Lisp implementation by Mario S. Mommer, Maxima interface by Andrej Vodopivec. Zip file retrieved from the old Maxima wiki many years ago. This commit includes the package as composed by Vodopivec verbatim, in addition to a basic Texinfo documentation file and regression test script. Also update build machinery to execute make in share/nelder_mead to run makeinfo and generate documentation index. --- configure.ac | 1 + share/Makefile.am | 2 +- share/nelder_mead/.gitignore | 6 + share/nelder_mead/COPYING.txt | 20 + share/nelder_mead/Makefile.am | 15 + share/nelder_mead/README.txt | 188 +++++++++ share/nelder_mead/defpackage.lisp | 5 + share/nelder_mead/la.lisp | 95 +++++ share/nelder_mead/load-neldermead.lisp | 9 + share/nelder_mead/nelder_mead.mac | 1 + share/nelder_mead/nelder_mead.system | 10 + share/nelder_mead/nelder_mead.texi | 90 ++++ share/nelder_mead/nelder_mead.txt | 22 + share/nelder_mead/neldermead.asd | 11 + share/nelder_mead/neldermead.lisp | 707 ++++++++++++++++++++++++++++++++ share/nelder_mead/nm-maxima.lisp | 12 + share/nelder_mead/rtest_nelder_mead.mac | 19 + src/share-subdirs.lisp | 2 +- 18 files changed, 1213 insertions(+), 2 deletions(-) create mode 100644 share/nelder_mead/.gitignore create mode 100644 share/nelder_mead/COPYING.txt create mode 100644 share/nelder_mead/Makefile.am create mode 100644 share/nelder_mead/README.txt create mode 100644 share/nelder_mead/defpackage.lisp create mode 100644 share/nelder_mead/la.lisp create mode 100644 share/nelder_mead/load-neldermead.lisp create mode 100644 share/nelder_mead/nelder_mead.mac create mode 100644 share/nelder_mead/nelder_mead.system create mode 100644 share/nelder_mead/nelder_mead.texi create mode 100644 share/nelder_mead/nelder_mead.txt create mode 100644 share/nelder_mead/neldermead.asd create mode 100644 share/nelder_mead/neldermead.lisp create mode 100644 share/nelder_mead/nm-maxima.lisp create mode 100644 share/nelder_mead/rtest_nelder_mead.mac diff --git a/configure.ac b/configure.ac index 262f6e498..9b0a8b5fc 100644 --- a/configure.ac +++ b/configure.ac @@ -1242,6 +1242,7 @@ plotting/mgnuplot share/Makefile demo/Makefile plotting/Makefile locale/Makefile share/contrib/Makefile share/contrib/integration/Makefile \ share/contrib/maxima-odesolve/Makefile \ share/contrib/symplectic_ode/Makefile \ +share/nelder_mead/Makefile \ share/draw/Makefile share/logic/Makefile doc/info/es/include-maxima.texi \ src/lisp]) AC_OUTPUT diff --git a/share/Makefile.am b/share/Makefile.am index 27d14b8fd..f56124228 100644 --- a/share/Makefile.am +++ b/share/Makefile.am @@ -19,4 +19,4 @@ $(sharefiles) EXTRA_DIST = $(genericdirDATA) -SUBDIRS = contrib logic draw +SUBDIRS = contrib logic draw nelder_mead diff --git a/share/nelder_mead/.gitignore b/share/nelder_mead/.gitignore new file mode 100644 index 000000000..7082e42b0 --- /dev/null +++ b/share/nelder_mead/.gitignore @@ -0,0 +1,6 @@ +*.info +nelder_mead-index.lisp +texinfo.tex +*.t2p +*.html +*.pdf diff --git a/share/nelder_mead/COPYING.txt b/share/nelder_mead/COPYING.txt new file mode 100644 index 000000000..802658975 --- /dev/null +++ b/share/nelder_mead/COPYING.txt @@ -0,0 +1,20 @@ +Copyright (c) 2006, Mario S. Mommer + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/share/nelder_mead/Makefile.am b/share/nelder_mead/Makefile.am new file mode 100644 index 000000000..4e6c16212 --- /dev/null +++ b/share/nelder_mead/Makefile.am @@ -0,0 +1,15 @@ +all-local: info + +info: nelder_mead.info nelder_mead-index.lisp + +nelder_mead.info: nelder_mead.texi + makeinfo --force $< + +nelder_mead-index.lisp: nelder_mead.info + perl $(top_srcdir)/doc/info/build_index.pl $< > $@ + +info_TEXINFOS = nelder_mead.texi + +AM_MAKEINFOHTMLFLAGS = --no-split + +EXTRA_DIST = nelder_mead-index.lisp diff --git a/share/nelder_mead/README.txt b/share/nelder_mead/README.txt new file mode 100644 index 000000000..a5b5eb182 --- /dev/null +++ b/share/nelder_mead/README.txt @@ -0,0 +1,188 @@ + Common-Lisp implementations of the + ================================== + + grid-restrained and traditional + =============================== + + Nelder-Mead algorithms + ====================== + + + +Introduction +------------ + + These common lisp sources contain two variants of the Nelder-Mead + algorithm. The original algorithm [1] and a provably convergent, + reliable variant by A. Bürmen et al [4], called the "Grid Restrained + Nelder Mead Algorithm" (GRNMA). + + It should be mentioned that other, provably convergent variant exist + [2,3], which aren't included here. The only reasons are lack of + time, and the fact that the implemented variant does not require a + simple descent condition, putting it closer to the original. + + Other than that, and based on the article [4], the performance of + these methods seems to be about equal in terms of number of function + evaluations. As a side effect of the additional reliability, both + tend to be a lot more efficient than the original algorithm even + when it does not fail. In particular when the number of dimensions + increases. As a test, one might try the GRNM, + + (grnm-optimize #'standard-quadratic + (make-array 30 :initial-element 1.0d0) :verbose t) + + and compare with the original Nelder-Mead, + + (nm-optimize #'standard-quadratic + (make-array 30 :initial-element 1.0d0) :verbose t) + + and observe the difference in number of function evaluations (the + last value returned). + + (This exercise also serves to illustrate the overall deficiencies of + direct search algorithms when applied to higher dimensional + problems.) + + This software is provided under the MIT licence; see COPYING.txt for + details. + +Usage +----- + + This implementation of the grid restrained Nelder-Mead algorithm + expects at least two parameters: the objective function, and an + initial guess. It returns four values; + + * the minimizer, + * the minimum, + * the last simplex, + * and the number of function evaluations. + + The objective function should accept a double-float array as its + argument. + + The initial guess will usually be an array of double-float numbers + (but instead an object of the class NM-SIMPLEX may also be provided; + see below). For example, + + (grnm-optimize #'rosenbrock #(40.0d0 40.0d0)) + + finds the minimum of the Rosenbrock function, and + + (grnm-optimize #'standard-quadratic + (make-array 30 :initial-element 1.0d0)) + + finds the minimum of the standard quadratic function in 30 + dimensions. + + A few keyword arguments can customize the behavior. These are + + :verbose (default: NIL) + :converged-p (default: (burmen-et-al-convergence-test)) + :max-function-calls (default: NIL; => as many as needed) + + + :verbose (default: NIL) + + Pass T here if you want to see some progress report. The amount of + output can be controlled by setting *verbose-level* to 1 or 2. The + difference is that with 1 (the default) only the best value of the + simplex is shown, while with 2 the whole simplex is printed on + each iteration. + + :converged-p (default: (burmen-et-al-convergence-test)) + + The burmen-et-al-convergence-test is, as the name suggests, the + convergence test used in the article by Burmen et al. It accepts a + few parameters: tol-x, tol-f and rel; please see the article for + further details. + + Another convergence criterion that can be given is (pp-volume-test + ), which returns true once the parallelogram(!) spanned by the + vertices of the simplex has a volume lower than to the power + of N, where N is the dimension of the problem. This is a rather + expensive test, as it involves computing a QR decomposition of an N + by N matrix on each call. + + If you are not in the mood of taking prisoners, you might as well + pass (constantly NIL) as the convergence criterion. This has as a + consequence that the iteration continues until the simplex + collapses, which in floating-point arithmetic happens in finite + time. The grid restrained Nelder-Mead algorithm should have + converged by then. + + :max-function-calls (default: NIL) + + Maximum number of objective function evaluations. After that many + function evaluations, this implementation of the algorithm will + declare convergence to have occurred. + + The actual number of function calls might be slightly larger (at + most by N), as the relevant condition is only checked in certain + situations. + +Utilities/Misc +-------------- + + NM-optimize + + Apart from the grid-restrained Nelder Mead algorithm, the + traditional variant is also provided. The corresponding function + is named NM-optimize, and its usage is the same as for + GRNM-optimize. + + The only difference is that :max-function-calls has a default + value of 100000. Otherwise the algorithm might well iterate + forever. + + Initial-simplex :displace + + Constructs an initial simplex with the double-float array as on of its corners. + + The displacement can be a an array of double floats or a + number. In the first case, the additional vertices of the simplex + are build by adding to each component of the the + corresponding component in . If it is a number, then + the additional vertices are build by adding to each + component of the . + +Tips & Tricks +------------- + + * One can try to solve constrained optimization problems by returning + + MOST-POSITIVE-DOUBLE-FLOAT + + whenever the objective function is called with an argument that + violates the constraints. Mathematically, this falls out of the + theory, but it works most of the time. + + * If your objective function is noisy or not smooth (for instance, if + its first derivatives are not continuous) it is a good idea to + restart the algorithm. Remember that convergence is only + guaranteed if the objective function is at least C^1. + +References +---------- + +[1] J.A. Nelder and R. Mead, "A simplex method for function + minimization," The Computer Journal, vol. 7, pp. 308-313, 1965. + + +[2] P. Tseng, "Fortified-descent simplicial search method: A general + approach," SIAM Journal on Optimization, vol. 10, pp. 269-288, + 1999. + +[3] C.J. Price, I.D. Coope, and D. Byatt, "A convergent variant of the + Nelder-Mead algorithm," Journal of Optimization Theory and + Applications, vol. 113, pp. 5-19, 2002. + +[4] A. Bürmen, J. Puhan and T. Tuma, "Grid Restrained Nelder-Mead + Algorithm", Computational Optimization and Applications, vol. + 34, no. 3, pp. 359 - 375, 2006. + + +-------- +Mario S. Mommer, November 2006 diff --git a/share/nelder_mead/defpackage.lisp b/share/nelder_mead/defpackage.lisp new file mode 100644 index 000000000..5bff3ab9c --- /dev/null +++ b/share/nelder_mead/defpackage.lisp @@ -0,0 +1,5 @@ +(defpackage #:neldermead + (:use :cl) + (:export :nm-optimize :grnm-optimize :initial-simplex + :rosenbrock :standard-quadratic :pp-volume-test + :burmen-et-al-convergence-test)) \ No newline at end of file diff --git a/share/nelder_mead/la.lisp b/share/nelder_mead/la.lisp new file mode 100644 index 000000000..63e91a661 --- /dev/null +++ b/share/nelder_mead/la.lisp @@ -0,0 +1,95 @@ +(in-package :neldermead) + +;; Simple matrix-vector facility. Optimized for staying out of the way. +(defmacro with-matrix-dimensions (varforms &body body) + (let ((form `(progn ,@body))) + (dolist (vform varforms) + (setf form + `(destructuring-bind ,(car vform) + (array-dimensions ,(cadr vform)) + (declare (fixnum ,@(car vform)) + (ignorable ,@(car vform))) + ,form))) + form)) + +(defun ip (v w &optional (start 0)) + (let ((l (length v)) + (ax 0.0d0)) + (dotimes (i (- l start)) + (incf ax (* (aref v (+ i start)) + (aref w (+ i start))))) + ax)) + +(defun norm (v &optional (start 0)) + (sqrt (ip v v start))) + +(defun make-matrix (m n) + (make-array (list m n) + :element-type 'double-float + :initial-element 0.0d0)) + +(defun make-vector (n &key (initial-element 0.0d0)) + (make-array n + :element-type 'double-float + :initial-element initial-element)) + +(defun v*c (vector constant) + (let* ((n (length vector)) + (res (make-vector n))) + (dotimes (i n) + (setf (aref res i) + (* constant (aref vector i)))) + res)) + +(defun v+w*c (v w c) + (assert (= (length v) (length w))) + (let* ((n (length v)) + (res (make-vector n))) + (dotimes (i n) + (setf (aref res i) + (+ (aref v i) + (* c (aref w i))))) + res)) + +(defun qr-factorization (mat &key (with-q t)) + (declare (type (simple-array double-float (* *)) mat) + (optimize (speed 3) (safety 0))) + (with-matrix-dimensions (((m n) mat)) + (let ((q (when with-q + (let ((pq (make-matrix m n))) + (declare (type (simple-array double-float (* *)) pq)) + (dotimes (i (min m n)) + (setf (aref pq i i) 1.0d0)) + pq)))) + + (loop for i from 0 below (- m 1) do + (loop for j from (+ i 1) below m do + (let ((a (aref mat i i)) + (b (aref mat j i))) + (when (/= b 0.0d0) + (let* ((r (sqrt (+ (* a a) (* b b)))) + (c (/ a r)) + (s (/ b r))) + + (loop for k from i below m do + (let ((olda (aref mat i k)) + (oldb (aref mat j k))) + + (setf (aref mat i k) + (+ (* c olda) (* s oldb)) + (aref mat j k) + (+ (* (- s) olda) (* c oldb))))) + + (when with-q + (let ((q q)) + (declare (type (simple-array double-float (* *)) q)) + (loop for k from 0 below m do + (let ((olda (aref q k i)) + (oldb (aref q k j))) + + (setf (aref q k i) + (+ (* c olda) (* s oldb)) + (aref q k j) + (+ (* (- s) olda) (* c oldb)))))))))))) + + (values mat q)))) diff --git a/share/nelder_mead/load-neldermead.lisp b/share/nelder_mead/load-neldermead.lisp new file mode 100644 index 000000000..228c63669 --- /dev/null +++ b/share/nelder_mead/load-neldermead.lisp @@ -0,0 +1,9 @@ +(in-package #-gcl #:maxima #+gcl "MAXIMA") + +#+ecl ($load "lisp-utils/defsystem.lisp") + +(load (merge-pathnames (make-pathname :name "nelder_mead" :type "system") + #-gcl *load-pathname* + #+gcl sys:*load-pathname*)) + +(mk:oos "nelder_mead" :compile) diff --git a/share/nelder_mead/nelder_mead.mac b/share/nelder_mead/nelder_mead.mac new file mode 100644 index 000000000..18eef930a --- /dev/null +++ b/share/nelder_mead/nelder_mead.mac @@ -0,0 +1 @@ +load("nelder_mead/load-neldermead.lisp"); diff --git a/share/nelder_mead/nelder_mead.system b/share/nelder_mead/nelder_mead.system new file mode 100644 index 000000000..fa0efb04f --- /dev/null +++ b/share/nelder_mead/nelder_mead.system @@ -0,0 +1,10 @@ + +(mk:defsystem nelder_mead + :source-pathname (maxima::maxima-load-pathname-directory) + :binary-pathname (maxima::maxima-objdir "share" "contrib" "nelder_mead") + :source-extension "lisp" + :components + ((:file "defpackage") + (:file "la.lisp") + (:file "neldermead.lisp") + (:file "nm-maxima.lisp"))) diff --git a/share/nelder_mead/nelder_mead.texi b/share/nelder_mead/nelder_mead.texi new file mode 100644 index 000000000..4164cf1e4 --- /dev/null +++ b/share/nelder_mead/nelder_mead.texi @@ -0,0 +1,90 @@ +\input texinfo + +@setfilename nelder_mead.info +@settitle Package nelder_mead + +@ifinfo +@macro var {expr} +<\expr\> +@end macro +@end ifinfo + +@dircategory Mathematics/Maxima +@direntry +* Package nelder_mead: (maxima)Maxima share package implementing the Nelder-Mead minimization algorithm +@end direntry + +@node Top, Introduction to package nelder_mead, (dir), (dir) +@top +@menu +* Introduction to package nelder_mead:: +* Definitions for package nelder_mead:: +* Function and variable index:: +@end menu +@chapter Package nelder_mead + +@node Introduction to package nelder_mead, Definitions for package nelder_mead, Top, Top +@section Introduction to package nelder_mead + +Package @code{nelder_mead} implements the Nelder-Mead minimization algorithm [1], +also known as the polytope or amoeba method. + +The Nelder-Mead algorithm is a derivative-free minimization algorithm; +only evaluations of the objective function are required. + +@code{nelder_mead} implements the "grid-restrained" Nelder-Mead algorithm published by A. Bürmen et al. [2], +and implemented in Common Lisp by Mario S. Mommer. +Thanks to Andrej Vodopivec for the Maxima interface to the Common Lisp code. + +References + +[1] J.A. Nelder and R. Mead, "A simplex method for function + minimization," The Computer Journal, vol. 7, pp. 308-313, 1965. + +[2] A. Bürmen, J. Puhan and T. Tuma, "Grid Restrained Nelder-Mead + Algorithm", Computational Optimization and Applications, vol. + 34, no. 3, pp. 359 - 375, 2006. + +@node Definitions for package nelder_mead, Function and variable index, Introduction to package nelder_mead, Top +@section Definitions for package nelder_mead + +@deffn {Function} nelder_mead (@var{obj}, @var{vars}, @var{init}) + +Returns an approximate minimum of the objective function @var{obj}, +as a function of the variables @var{vars}, +starting at the initial point @var{init}. + +The objective function may be discontinuous, +but if it is continuous, +the algorithm ("grid-restrained" Nelder-Mead) +is provably convergent. +The objective function need not be differentiable; +derivatives are not computed, not even approximately. + +Examples: + +@c ===beg=== +@c load ("nelder_mead") $ +@c nelder_mead (if x<0 then -x else x^2, [x], [4]); +@c f(x) := if x < 0 then -x else x^2 $ +@c nelder_mead (f, [x], [4]); +@c nelder_mead (x^4 + y^4 - 2*x*y - 4*x - 3*y, [x, y], [2, 2]); +@c ===end=== +@example +(%i1) load ("nelder_mead") $ +(%i2) nelder_mead (if x<0 then -x else x^2, [x], [4]); +(%o2) [x = 9.536387892694628E-11] +(%i3) f(x) := if x < 0 then -x else x^2 $ +(%i4) nelder_mead (f, [x], [4]); +(%o4) [x = 9.536387892694628E-11] +(%i5) nelder_mead (x^4 + y^4 - 2*x*y - 4*x - 3*y, [x, y], [2, 2]); +(%o5) [x = 1.157212489168102, y = 1.099342680267472] +@end example +@end deffn + +@node Function and variable index, , Definitions for package nelder_mead, Top +@appendix Function and variable index +@printindex fn +@printindex vr + +@bye diff --git a/share/nelder_mead/nelder_mead.txt b/share/nelder_mead/nelder_mead.txt new file mode 100644 index 000000000..3ea7f2d58 --- /dev/null +++ b/share/nelder_mead/nelder_mead.txt @@ -0,0 +1,22 @@ + + +Nelder-Mead minimization method + +Postby andrejv » Tue Jun 16, 2009 7:03 am +This package contains the Nelder-Mead minimization method. + +The method was written in lisp by Mario S. Mommer (http://prxq.wordpress.com/2006/11/05/grid-restrained-nelder-mead/) and is available under the MIT license. + +This package requires Maxima version 5.18.2. + + (%i1) load("nelder_mead/nelder_mead.mac")$ + (%i2) nelder_mead(if x<0 then -x else x^2, [x], [4]); + (%o2) [x = 9.536387892694628e-11] + (%i3) f(x) := if x<0 then -x else x^2$ + (%i4) nelder_mead(f, [x], [4]); + (%o4) [x = 9.536387892694628e-11] + (%i5) nelder_mead(f(x), [x], [4]); + (%o5) [x = 9.536387892694628e-11] + (%i6) nelder_mead(x^4+y^4-2*x*y-4*x-3*y, [x,y], [2,2]); + (%o6) [x = 1.157212489168102, y = 1.099342680267472] + diff --git a/share/nelder_mead/neldermead.asd b/share/nelder_mead/neldermead.asd new file mode 100644 index 000000000..b86ab7437 --- /dev/null +++ b/share/nelder_mead/neldermead.asd @@ -0,0 +1,11 @@ +; -*- lisp -*- +(defpackage #:neldermead-asd + (:use :cl :asdf)) + +(in-package :neldermead-asd) + +(defsystem :neldermead + :serial t + :components ((:file "defpackage") + (:file "la") + (:file "neldermead"))) \ No newline at end of file diff --git a/share/nelder_mead/neldermead.lisp b/share/nelder_mead/neldermead.lisp new file mode 100644 index 000000000..94199eb5d --- /dev/null +++ b/share/nelder_mead/neldermead.lisp @@ -0,0 +1,707 @@ +(in-package :neldermead) + +(defclass cached-simplex-data () + ((pseudopivot :initform nil) + (q-factor :initform nil) + (side-vectors :initform nil) + (r-factor :initform nil) + (dv :initform nil + :initarg :dv + :accessor dv) + (dmin :initform nil + :initarg :dmin + :accessor dmin) + (best-before-reshape :initform nil + :initarg :best-before-reshape + :accessor best-before-reshape))) + +(defclass nm-simplex () + ((x :accessor x + :initarg :x) + (fx :accessor fx + :initarg :fx) + (pmap :accessor pmap + :initarg :pmap) + (data :accessor data + :initform nil))) + +(defclass grid () + ((z :accessor grid-z :initarg :z) + (delta :accessor delta :initarg :delta))) + +(defvar *verbose-level* 1) + +(defmethod print-object ((s nm-simplex) stream) + (format stream "#")) + +(defmethod xk ((s nm-simplex) k) + (aref (x s) (aref (pmap s) k))) + +(defmethod (setf xk) (nv (s nm-simplex) k) + (setf (aref (x s) (aref (pmap s) k)) nv)) + +(defmethod fk ((s nm-simplex) k) + (aref (fx s) (aref (pmap s) k))) + +(defmethod (setf fk) (nv (s nm-simplex) k) + (setf (aref (fx s) (aref (pmap s) k)) nv)) + +(defmethod sort-simplex ((s nm-simplex)) + (sort (pmap s) #'< :key #'(lambda (k) (aref (fx s) k))) + s) + +(defmethod dimension ((s nm-simplex)) + (- (length (pmap s)) 1)) + +;; simple additive simplex generator +(defun initial-simplex (x0 + &key (displace 0.1d0)) + (let* ((n (length x0)) + (x (make-array (+ n 1))) + (pmap (make-array (+ n 1) :element-type 'fixnum))) + + (setf (aref x 0) x0 + (aref pmap 0) 0) + + (dotimes (k n) + (let ((xk (copy-seq x0))) + + (incf (aref xk k) + (if (numberp displace) displace + (aref displace k))) + (setf (aref x (+ k 1)) xk + (aref pmap (+ k 1)) (+ k 1)))) + + (make-instance 'nm-simplex + :x x :fx nil :pmap pmap))) + +;; The simplex generator from the article, more or less. +(defun default-initial-simplex (x0) + (initial-simplex x0 + :displace + (map 'vector + #'(lambda (v) + (max 0.00025d0 + (abs (* v 0.05d0)))) + x0))) + +(defmethod maybe-fill-simplex ((s nm-simplex) f) + (if (fx s) s + (let* ((n (length (x s))) + (fx (make-array n :element-type 'double-float))) + (loop :for i :from 0 :below n :do + (setf (aref fx i) (funcall f (xk s i)))) + + (setf (fx s) fx) + s))) + +;; substitutes the worst point with x/fx. Assumes simplex is sorted. +(defmethod improve ((s nm-simplex) x fx) + (let ((last (length x))) + (setf (xk s last) x + (fk s last) fx + + (data s) nil) + + (sort-simplex s))) + +(defmethod cached-slot ((s nm-simplex) slot computer) + (if (and (data s) (slot-value (data s) slot)) + (slot-value (data s) slot) + (setf (data s) (if (data s) (data s) + (make-instance 'cached-simplex-data)) + (slot-value (data s) slot) (funcall computer)))) + +(defmethod pseudopivot ((s nm-simplex)) + (cached-slot + s 'pseudopivot + + #'(lambda () + (let* ((n (- (length (pmap s)) 1)) + (xbar (make-array n + :element-type 'double-float + :initial-element 0.0d0))) + + (dotimes (i n) + (setf xbar (v+w*c xbar (xk s (+ i 1)) (/ 1.0d0 n)))) + + xbar)))) + +(defmethod side-vectors ((s nm-simplex)) + (cached-slot + s 'side-vectors + + #'(lambda () + (let* ((n (- (length (pmap s)) 1)) + (sv (make-array n))) + + (dotimes (i n) + (setf (aref sv i) (v+w*c (xk s (+ i 1)) (xk s 0) -1))) + + sv)))) + +(defun simplex-qr-thing (sidev) + (let* ((n (length sidev)) + (norms (make-array n :element-type 'double-float)) + (pmap (make-array n :element-type 'fixnum)) + (mat (make-array (list n n) :element-type 'double-float))) + + (dotimes (i n) + (setf (aref pmap i) i + (aref norms i) (norm (aref sidev i)))) + + (sort pmap #'> :key #'(lambda (k) (aref norms k))) + + (dotimes (i n) + (dotimes (j n) + (setf (aref mat i j) + (aref (aref sidev (aref pmap j)) i)))) + + (multiple-value-bind (r q) (qr-factorization mat) + (values q r pmap)))) + +(defun qrthing-closure (s n) + #'(lambda () + (multiple-value-bind (q r p) + (simplex-qr-thing (side-vectors s)) + (setf (slot-value (data s) 'q-factor) q + (slot-value (data s) 'r-factor) r) + + (elt (list q r p) n)))) + +(defmethod q-factor ((s nm-simplex)) + (cached-slot s 'q-factor + + (qrthing-closure s 0))) + +(defmethod r-factor ((s nm-simplex)) + (cached-slot s 'r-factor + (qrthing-closure s 1))) + +;; A Nelder-Mead iteration. +(defun nm-iteration + (simplex f &key + verbose + (gamma_reflect 1.0d0) + (gamma_expand 2.0d0) + (gamma_outer_contraction 0.5d0) + (gamma_inner_contraction -0.5d0) + (gamma_shrink 0.5d0)) + + + (let* ((n (- (length (x simplex)) 1)) + + (x_cb (make-array n + :element-type 'double-float + :initial-element 0.0d0)) + (x_cb-x_n (make-array n + :element-type 'double-float + :initial-element 0.0d0))) + + (labels ((newpoint (gamma) + (let ((np (v+w*c x_cb x_cb-x_n gamma))) + (values np (funcall f np)))) + + (accept (xx ff) + (improve simplex xx ff)) + + (shrink () + (let ((x0 (xk simplex 0))) + (loop for i from 1 to n do + (let* ((newx (v+w*c (v*c x0 + (- 1.0d0 gamma_shrink)) + (xk simplex i) gamma_shrink)) + (newf (funcall f newx))) + + (setf (xk simplex i) newx + (fk simplex i) newf) + + (sort-simplex simplex)))))) + + ;; compute centroid + (dotimes (i n) + (setf x_cb (v+w*c x_cb (xk simplex i) (/ 1.0d0 n)))) + + (setf x_cb-x_n (v+w*c x_cb (xk simplex n) -1.0d0)) + + ;; 2. Reflect + (multiple-value-bind (xr fr) + (newpoint gamma_reflect) + (if (and (<= (fk simplex 0) fr) (< fr (fk simplex (- n 1)))) + (accept xr fr) + ;; 3. expand + (if (< fr (fk simplex 0)) + (multiple-value-bind (xe fe) + (newpoint gamma_expand) + (if (< fe fr) + (accept xe fe) + (accept xr fr))) + ;; 4. contract or shrink + (if (<= (fk simplex (- n 1)) fr) + (if (< fr (fk simplex n)) + ;;outer contraction + (multiple-value-bind (xc fc) + (newpoint gamma_outer_contraction) + (if (or (= n 2) + (<= fc fr)) ;; 1D shrink is + ;; equivalent to + ;; acceptance of this + ;; point + (accept xc fc) + (shrink))) + ;; inner contraction + (multiple-value-bind (xcc fcc) + (newpoint gamma_inner_contraction) + (if (< fcc (fk simplex n)) (accept xcc fcc) + (shrink)))) + + (shrink))))) + + (when verbose (format t "~S~%" simplex)) + simplex))) + +;; The test returns true if the volume of the paralelepiped spanned by +;; the vertices of the simplex is smaller than that of an n-cube of +;; side cside. +(defun pp-volume-test (cside) + #'(lambda (simplex) + (let* ((mat (r-factor simplex))) + + (let ((det 1.0d0)) + (dotimes (i (dimension simplex)) + (setf det (* det (aref mat i i)))) + + (< (abs det) (expt cside (dimension simplex))))))) + +(defun nm-optimize (objective-function initial-guess &key + (max-function-calls 100000) + (convergence-p (burmen-et-al-convergence-test)) + verbose) + + (let ((simplex (if (typep initial-guess 'nm-simplex) + initial-guess + (default-initial-simplex initial-guess))) + (fvcount 0)) + (labels ((rigged-f (v) + (incf fvcount) + (funcall objective-function v)) + (converged-p (s) + (or (funcall convergence-p s) + (> fvcount max-function-calls)))) + + (when verbose + (format t "Initial simplex: ~%~A~%---~%" simplex)) + + (maybe-fill-simplex simplex #'rigged-f) + + (loop :until (converged-p simplex) + :do + (nm-iteration simplex #'rigged-f :verbose verbose)) + + (values (xk simplex 0) (fk simplex 0) simplex fvcount)))) + +(defmethod restrict ((grid grid) point) + (let ((new (copy-seq point)) + (n (length point)) + (delta (delta grid)) + (z (grid-z grid))) + + (dotimes (i n) + (setf (aref new i) + (+ (* (aref delta i) + (floor + (+ (/ (- (aref new i) + (aref z i)) + (aref delta i)) + 0.5d0))) + (aref z i)))) + + new)) + + +;; Some parameters. Look at Burmen et al for further details. +(defparameter *psi* 1.0d-6) +(defparameter *biglambda* (/ 0.5d0 double-float-epsilon)) +(defparameter *tau-r* (* 2.0d0 double-float-epsilon)) +(defparameter *tau-a* (expt least-positive-double-float (/ 1.0d0 3.0d0))) +(defparameter *smalllambda* 2) + +(defvar *breakdown*) + +(defmethod maybe-reshape ((s nm-simplex) (g grid) ff &key force) + (let ((n (length (side-vectors s))) + (biglambda *biglambda*) + (smalllambda *smalllambda*) + (reshaped-p nil)) + + (labels ((degenerate-p () + (let* ((r (r-factor s)) + (ax (loop for i from 0 below n + minimizing (abs (aref r i i))))) + + (< ax + (/ (* *psi* + (norm (delta g)) + (sqrt (float n 1.0d0))) + 2))))) + + (when (or force (degenerate-p)) + + (setf reshaped-p t) + + (let ((r (r-factor s)) + (q (q-factor s)) + (|det| 1.0d0) + (dv (make-array (+ n 1))) + (dmin-norm nil) + (dmin nil) + (|Delta| (norm (delta g)))) + + (setf (aref dv 0) (xk s 0)) + + (dotimes (i n) + (let* ((di (make-array n :element-type 'double-float)) + (rii (aref r i i)) + (sgn[rii] (if (>= rii 0) 1 -1)) + (|rii| (abs rii)) + (quot (* (sqrt (float n 0.0d0)) + |Delta| + 0.5d0)) + (minc (min |rii| + (* biglambda quot))) + (maxc (max (* smalllambda quot) minc)) + (dfkt (* sgn[rii] maxc))) + + (setf |det| (* |det| |rii|)) + + (dotimes (j n) + (setf (aref di j) + (* dfkt (aref q j i)))) + + (when (or (not dmin-norm) (< (abs dfkt) dmin-norm)) + (setf dmin-norm (abs dfkt) + dmin di)) + + (setf (aref dv (+ i 1)) di) + + (let* ((nxi (restrict g (v+w*c (xk s 0) di 1.0d0))) + (fxi (funcall ff nxi))) + + (setf (xk s (+ i 1)) nxi + (fk s (+ i 1)) fxi)))) + + ;; Looks like a very extreme situation, but can actually + ;; happen. + (when (= |det| 0.0d0) + (setf *breakdown* t)) + + ;; Most cached data is invalid, but better keep the reshape + ;; data, which might be needed for shrinking the simplex. + (setf (data s) + (make-instance 'cached-simplex-data + :dv dv + :dmin dmin + :best-before-reshape (fk s 0))) + + (sort-simplex s)))) + + reshaped-p)) + +;;; Nelder Mead iteration variant from Burmen et al. Does not shrink, +;;; "failing" instead. Acceptance criteria for the contraction points +;;; are stricter. +;;; +;;; Iterates are restricted to a grid. +(defun nm-iteration-burmen-et-al + (simplex f grid &key + verbose + (gamma_reflect 1.0d0) + (gamma_expand 2.0d0) + (gamma_outer_contraction 0.5d0) + (gamma_inner_contraction -0.5d0)) + + + (let ((n (- (length (x simplex)) 1)) + (failure t)) + + (let ((x_cb (make-array n + :element-type 'double-float + :initial-element 0.0d0)) + (x_cb-x_n (make-array n + :element-type 'double-float + :initial-element 0.0d0))) + + (labels ((newpoint (gamma) + (let ((np (restrict grid + (v+w*c x_cb x_cb-x_n gamma)))) + (values np (funcall f np)))) + + (accept (xx ff) + (improve simplex xx ff) + (setf failure nil))) + + ;; compute centroid + (dotimes (i n) + (setf x_cb (v+w*c x_cb (xk simplex i) (/ 1.0d0 n)))) + + (setf x_cb-x_n (v+w*c x_cb (xk simplex n) -1.0d0)) + + ;; 2. Reflect + (multiple-value-bind (xr fr) + (newpoint gamma_reflect) + (if (and (<= (fk simplex 0) fr) (< fr (fk simplex (- n 1)))) + (accept xr fr) + ;; 3. expand + (if (< fr (fk simplex 0)) + (multiple-value-bind (xe fe) + (newpoint gamma_expand) + (if (< fe fr) + (accept xe fe) + (accept xr fr))) + ;; 4. contract or fail ABurmen & al have swapped + ;; inner and outers - maybe a typo. (?) No, just a + ;; different variant. + (if (<= (fk simplex (- n 1)) fr) + (if (< fr (fk simplex n)) + ;;outer contraction + (multiple-value-bind (xc fc) + (newpoint gamma_outer_contraction) + (if (<= fc (fk simplex (- n 1))) + (accept xc fc) + )) + ;; inner contraction + (multiple-value-bind (xcc fcc) + (newpoint gamma_inner_contraction) + (if (< fcc (fk simplex (- n 1))) + (accept xcc fcc)))))))))) + + (sort-simplex simplex) + (when verbose (format t "~S~%" simplex)) + (values simplex failure))) + +;; Convergence test used in the article +(defun burmen-et-al-convergence-test (&key + (tol-x 1.0d-8) + (tol-f 1.0d-15) + (rel 1.0d-15)) + + #'(lambda (s) + (let* ((sv (side-vectors s)) + (fdiff (abs (- (fk s 0) (fk s (dimension s))))) + (vijmax (loop for v across sv maximizing + (loop for x across v maximizing + (abs x)))) + (xbest (xk s 0)) + (|xi|max (loop for z across xbest maximizing (abs z)))) + + (and (< fdiff (max tol-f (* rel (fk s 0)))) + (< vijmax (max tol-x (* rel |xi|max))))))) + + + +(defun regrid (g x1 dmin fkt) + (let* ((dmin (v*c dmin fkt)) + (n (length dmin)) + (|dmin| (norm dmin)) + (delta (delta g)) + (newd (copy-seq delta))) + + (setf (grid-z g) x1) + + (dotimes (i (length dmin)) + (setf (aref newd i) + + (max (min (max (/ (abs (aref dmin i)) + (* *smalllambda* 250 n)) + + (/ |dmin| + (* *smalllambda* 250 (expt n (/ 3.0d0 2.0d0))))) + (aref delta i)) + (* *tau-r* (aref x1 i)) + + *tau-a*))) + + (setf (delta g) newd))) + +(defun deep-shrink (f s g gamma_s convergence-p verbose) + (let ((bbr (fk s 0))) + (unless (dv (data s)) + (maybe-reshape s g f :force t)) + + (let* ((dv (dv (data s))) + (n (dimension s)) + (dmin (dmin (data s))) + (|dmin| (norm dmin)) + (converged-p nil) + (l 1)) ;; This does the same as in the article + + (loop :until (or (< (fk s 0) bbr) + (setf converged-p (funcall convergence-p s))) + :do (let ((fkt (* (expt gamma_s (floor l 2)) + (expt -1 l)))) + + (when verbose + (format t "Shrinking by factor: ~A~%" fkt)) + + (when (= 0.0d0 fkt) + ;; You've got no simplex anymore + (setf *breakdown* t)) + + (when (< (* (abs fkt) |dmin|) + (* (/ *smalllambda* 2) + (sqrt n) + (norm (delta g)))) + + (regrid g (xk s 0) dmin fkt)) + + (setf (xk s 0) (aref dv 0) + (fk s 0) bbr) + + (loop for k from 1 below (length (x s)) do + (setf (xk s k) (restrict g + (v+w*c (aref dv 0) + (aref dv k) fkt)) + (fk s k) (funcall f (xk s k)))) + + (sort-simplex s) + (incf l))) + + (setf (data s) nil) + + converged-p))) + +(defun grnm-optimize (objective-function initial-guess &key + max-function-calls + (convergence-p (burmen-et-al-convergence-test)) + verbose) + + (let ((simplex (if (typep initial-guess 'nm-simplex) + initial-guess + (default-initial-simplex initial-guess))) + (fvcount 0) + (gamma_s 0.5d0) + (*breakdown* nil)) + + (labels ((rigged-f (v) + (incf fvcount) + (funcall objective-function v)) + (converged-p (s) + (or (funcall convergence-p s) + *breakdown* + (and max-function-calls (> fvcount max-function-calls))))) + (when verbose + (format t "Initial simplex: ~%~A~%---~%" simplex)) + + (maybe-fill-simplex simplex #'rigged-f) + + (prog (failure xbest fbest pbest reshaped-p + + (grid (make-instance 'grid + :z (xk simplex 0) + :delta + (make-array (dimension simplex) + :initial-element + (/ (loop :for i :from 1 :to + (dimension simplex) + :minimizing + (norm + (v+w*c + (xk simplex 0) + (xk simplex i) -1.0d0))) + 10.0d0))))) + + + iterate ;; 1. + (setf (values simplex failure) + (nm-iteration-burmen-et-al simplex #'rigged-f grid + :verbose verbose)) + + ;; Small variation. We test for convergence here too. Depending on + ;; the convergence criterion, this might be spurious, so we have to + ;; reshape, etc. + (unless (or failure (converged-p simplex)) + (go iterate)) + + ;; 2. + (setf xbest (xk simplex 0) + fbest (fk simplex 0) + pbest (aref (pmap simplex) 0) + + reshaped-p (maybe-reshape simplex grid #'rigged-f)) + + + ;; 4. Here we look at the pseudo expand point + (let* ((pep (pseudopivot simplex)) + (xx (v+w*c pep + (v+w*c pep (xk simplex 0) -1.0d0) + 1.0d0)) ;; To do: clear this up + (fpep (rigged-f xx))) + + (if (<= fbest (min (fk simplex 0) fpep)) + (go deep-shrink) + (when (< fpep fbest) + ;; There is something subtle here. The pseudo-expand point + ;; is supposed to substitute the old (xk 0), which might + ;; have changed *position* during reshape, because of the + ;; simplex being sorted. + ;; + ;; So it is not + ;;(improve simplex xx fpep) + ;; here, but instead + + ;; Step 5. + (setf (aref (x simplex) pbest) xx + (aref (fx simplex) pbest) fpep + (data simplex) nil) + + (sort-simplex simplex) + + ;; Step 6. + (go iterate)))) + + deep-shrink + + (unless + (deep-shrink #'rigged-f simplex grid + + ;;(+ 0.1d0 (random 0.8d0)) + gamma_s + + #'converged-p verbose) + (go iterate)) + + end) + (values (xk simplex 0) (fk simplex 0) simplex fvcount)))) + + +;; Some sample functions to play around. + +(defun standard-quadratic (v) + (loop for i from 0 below (length v) + summing (expt (aref v i) 2))) + +(defun rosenbrock (v) + (let ((x (aref v 0)) + (y (aref v 1))) + + (+ (expt (- 1 x) 2) + (* 100 (expt (- y (expt x 2)) 2))))) diff --git a/share/nelder_mead/nm-maxima.lisp b/share/nelder_mead/nm-maxima.lisp new file mode 100644 index 000000000..ef9c0d22f --- /dev/null +++ b/share/nelder_mead/nm-maxima.lisp @@ -0,0 +1,12 @@ +(defun $nelder_mead (expr vars init) + (let* ((fun (coerce-float-fun expr vars)) + (fun1 (lambda (arr) + (mfuncall '$apply fun `((mlist simp) ,@(loop for i across arr collect i))))) + (init (make-array ($length init) :initial-contents (cdr ($float init))))) + (multiple-value-bind + (xk fk fv) (neldermead:grnm-optimize fun1 init :verbose nil) + (declare (ignore fk fv)) + `((mlist simp) ,@(mapcar #'(lambda (x y) `((mequal simp) ,x ,y)) + (cdr vars) + (loop for i across xk collect i)))))) + diff --git a/share/nelder_mead/rtest_nelder_mead.mac b/share/nelder_mead/rtest_nelder_mead.mac new file mode 100644 index 000000000..8da07fa2f --- /dev/null +++ b/share/nelder_mead/rtest_nelder_mead.mac @@ -0,0 +1,19 @@ +(if errcatch (fundef (nelder_mead)) = [] + then load ("nelder_mead"), + 0); +0; + +nelder_mead(if x<0 then -x else x^2, [x], [4]); +[x = 9.536387892694628e-11]; + +(f(x) := if x<0 then -x else x^2, 0); +0; + +nelder_mead(f, [x], [4]); +[x = 9.536387892694628e-11]; + +nelder_mead(f(x), [x], [4]); +[x = 9.536387892694628e-11]; + +nelder_mead(x^4+y^4-2*x*y-4*x-3*y, [x,y], [2,2]); +[x = 1.157212489168102, y = 1.099342680267472]; diff --git a/src/share-subdirs.lisp b/src/share-subdirs.lisp index 5c474dd2b..f57526176 100644 --- a/src/share-subdirs.lisp +++ b/src/share-subdirs.lisp @@ -19,5 +19,5 @@ ;; DO NOT EDIT THIS LIST. It is automatically ;; generated by configure. '( -"affine" "algebra" "algebra/charsets" "algebra/solver" "amatrix" "bernstein" "calculus" "cobyla" "cobyla/ex" "cobyla/lisp" "colnew" "colnew/ex1" "colnew/ex2" "colnew/ex3" "colnew/ex4" "colnew/lisp" "combinatorics" "contrib" "contrib/Eulix" "contrib/Grobner" "contrib/Zeilberger" "contrib/alt-display" "contrib/altsimp" "contrib/binsplit" "contrib/bitwise" "contrib/boolsimp" "contrib/coma" "contrib/diffequations" "contrib/diffequations/tests" "contrib/elliptic_curves" "contrib/elliptic_curves/figures" "contrib/format" "contrib/fresnel" "contrib/gentran" "contrib/gentran/man" "contrib/gentran/test" "contrib/gf" "contrib/integration" "contrib/levin" "contrib/lurkmathml" "contrib/maxima-odesolve" "contrib/maximaMathML" "contrib/mcclim" "contrib/namespaces" "contrib/noninteractive" "contrib/odes" "contrib/operatingsystem" "contrib/prim" "contrib/rand" "contrib/rkf45" "contrib/sarag" "contrib/smath" "contrib/state" "contrib/symplectic_ode" "contrib/trigtools" "contrib/unicodedata" "contrib/unit" "contrib/vector3d" "descriptive" "diff_form" "diff_form/tests" "diffequations" "distrib" "draw" "dynamics" "ezunits" "fftpack5" "fftpack5/lisp" "finance" "fourier_elim" "fractals" "graphs" "hompack" "hompack/lisp" "hypergeometric" "integequations" "integer_sequence" "integration" "lapack" "lapack/blas" "lapack/lapack" "lbfgs" "linearalgebra" "logic" "lsquares" "macro" "matrix" "minpack" "minpack/lisp" "misc" "mnewton" "multiadditive" "numeric" "numericalio" "odepack" "odepack/src" "orthopoly" "pdiff" "physics" "pslq" "pytranslate" "quantum" "simplex" "simplex/Tests" "simplification" "solve_rat_ineq" "solve_rec" "sound" "stats" "stringproc" "sym" "tensor" "test_batch_encodings" "to_poly_solve" "translators" "translators/m2mj" "trigonometry" "utils" "vector" "z_transform" +"affine" "algebra" "algebra/charsets" "algebra/solver" "amatrix" "bernstein" "calculus" "cobyla" "cobyla/ex" "cobyla/lisp" "colnew" "colnew/ex1" "colnew/ex2" "colnew/ex3" "colnew/ex4" "colnew/lisp" "combinatorics" "contrib" "contrib/Eulix" "contrib/Grobner" "contrib/Zeilberger" "contrib/alt-display" "contrib/altsimp" "contrib/binsplit" "contrib/bitwise" "contrib/boolsimp" "contrib/coma" "contrib/diffequations" "contrib/diffequations/tests" "contrib/elliptic_curves" "contrib/elliptic_curves/figures" "contrib/format" "contrib/fresnel" "contrib/gentran" "contrib/gentran/man" "contrib/gentran/test" "contrib/gf" "contrib/integration" "contrib/levin" "contrib/lurkmathml" "contrib/maxima-odesolve" "contrib/maximaMathML" "contrib/mcclim" "contrib/namespaces" "contrib/noninteractive" "contrib/odes" "contrib/operatingsystem" "contrib/prim" "contrib/rand" "contrib/rkf45" "contrib/sarag" "contrib/smath" "contrib/state" "contrib/symplectic_ode" "contrib/trigtools" "contrib/unicodedata" "contrib/unit" "contrib/vector3d" "descriptive" "diff_form" "diff_form/tests" "diffequations" "distrib" "draw" "dynamics" "ezunits" "fftpack5" "fftpack5/lisp" "finance" "fourier_elim" "fractals" "graphs" "hompack" "hompack/lisp" "hypergeometric" "integequations" "integer_sequence" "integration" "lapack" "lapack/blas" "lapack/lapack" "lbfgs" "linearalgebra" "logic" "lsquares" "macro" "matrix" "minpack" "minpack/lisp" "misc" "mnewton" "multiadditive" "nelder_mead" "numeric" "numericalio" "odepack" "odepack/src" "orthopoly" "pdiff" "physics" "pslq" "pytranslate" "quantum" "simplex" "simplex/Tests" "simplification" "solve_rat_ineq" "solve_rec" "sound" "stats" "stringproc" "sym" "tensor" "test_batch_encodings" "to_poly_solve" "translators" "translators/m2mj" "trigonometry" "utils" "vector" "z_transform" ))) -- 2.11.4.GIT