1
; libctl: flexible Guile-based control files for scientific software
2
; Copyright (C) 1998, 1999, 2000, 2001, 2002, Steven G. Johnson
4
; This library is free software; you can redistribute it and/or
5
; modify it under the terms of the GNU Lesser General Public
6
; License as published by the Free Software Foundation; either
7
; version 2 of the License, or (at your option) any later version.
9
; This library is distributed in the hope that it will be useful,
10
; but WITHOUT ANY WARRANTY; without even the implied warranty of
11
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12
; Lesser General Public License for more details.
14
; You should have received a copy of the GNU Lesser General Public
15
; License along with this library; if not, write to the
16
; Free Software Foundation, Inc., 59 Temple Place - Suite 330,
17
; Boston, MA 02111-1307, USA.
19
; Steven G. Johnson can be contacted at stevenj@alum.mit.edu.
21
; ****************************************************************
22
; The Nelder-Mead simplex algorithm for multidimensional minimization.
23
; See the simplex-minimize function, below.
25
(define (ax+by a x b y)
26
(define (cdr-null x) (if (null? x) '() (cdr x)))
27
(if (and (null? x) (null? y))
29
(let ((ax (if (null? x) 0 (* a (car x))))
30
(by (if (null? y) 0 (* b (car y)))))
31
(cons (+ ax by) (ax+by a (cdr-null x) b (cdr-null y))))))
33
(define (simplex-point x val) (cons x val))
34
(define simplex-point-x car)
35
(define simplex-point-val cdr)
37
(define (simplex-high s)
38
(car (sort s (lambda (s1 s2) (> (simplex-point-val s1)
39
(simplex-point-val s2))))))
40
(define (simplex-high2 s)
41
(cadr (sort s (lambda (s1 s2) (> (simplex-point-val s1)
42
(simplex-point-val s2))))))
43
(define (simplex-low s)
44
(car (sort s (lambda (s1 s2) (< (simplex-point-val s1)
45
(simplex-point-val s2))))))
47
(define (simplex-replace s s-old s-new)
50
(if (eq? (car s) s-old)
52
(cons (car s) (simplex-replace (cdr s) s-old s-new)))))
54
(define (simplex-sum-x s)
57
(ax+by 1 (simplex-point-x (car s)) 1 (simplex-sum-x (cdr s)))))
59
(define (simplex-centroid-x s)
60
(let ((sum (ax+by 1 (simplex-sum-x s)
61
-1 (simplex-point-x (simplex-high s)))))
62
(ax+by (/ (- (length s) 1)) sum 0.0 '())))
64
(define (simplex-shrink s-min f s)
67
(if (eq? s-min (car s))
68
(cons (car s) (simplex-shrink s-min f (cdr s)))
69
(let ((x (ax+by 0.5 (simplex-point-x s-min)
70
0.5 (simplex-point-x (car s)))))
71
(cons (simplex-point x (apply f x))
72
(simplex-shrink s-min f (cdr s)))))))
74
(define simplex-reflect-ratio 1.0)
75
(define simplex-expand-ratio 2.0)
76
(define simplex-contract-ratio 0.5)
78
(define (simplex-contract f s)
79
(let ((s-h (simplex-high s))
81
(x0 (simplex-centroid-x s)))
82
(let ((xc (ax+by (- 1 simplex-contract-ratio) x0
83
simplex-contract-ratio (simplex-point-x s-h))))
84
(let ((vc (apply f xc)))
85
(if (< vc (simplex-point-val s-h))
86
(simplex-replace s s-h (simplex-point xc vc))
87
(simplex-shrink s-l f s))))))
89
(define (simplex-iter f s)
90
(let ((s-h (simplex-high s))
91
(s-h2 (simplex-high2 s))
93
(x0 (simplex-centroid-x s)))
94
(let ((xr (ax+by (+ 1 simplex-reflect-ratio) x0
95
(- simplex-reflect-ratio) (simplex-point-x s-h))))
96
(let ((vr (apply f xr)))
97
(if (and (<= vr (simplex-point-val s-h2))
98
(>= vr (simplex-point-val s-l)))
100
(simplex-replace s s-h (simplex-point xr vr))
101
(if (< vr (simplex-point-val s-l))
102
(let ((xe (ax+by (- 1 simplex-expand-ratio) x0
103
simplex-expand-ratio xr)))
104
(let ((ve (apply f xe)))
106
(simplex-replace s s-h (simplex-point xr vr))
107
(simplex-replace s s-h (simplex-point xe ve)))))
108
(if (and (< vr (simplex-point-val s-h))
109
(> vr (simplex-point-val s-h2)))
110
(simplex-contract f (simplex-replace
111
s s-h (simplex-point xr vr)))
112
(simplex-contract f s))))))))
114
(define (simplex-iterate f s tol)
115
(let ((s-h (simplex-high s))
116
(s-l (simplex-low s)))
117
(if (<= (magnitude (- (simplex-point-val s-h) (simplex-point-val s-l)))
118
(* 0.5 tol (+ tol (magnitude (simplex-point-val s-h))
119
(magnitude (simplex-point-val s-l)))))
122
(print "extremization: best so far is " s-l "\n")
123
(simplex-iterate f (simplex-iter f s) tol)))))
125
(define (simplex-shift-x x i)
126
(let ((xv (list->vector x)))
127
(let ((xv-i (vector-ref xv i)))
128
(if (< (magnitude xv-i) 1e-6)
129
(vector-set! xv i 0.1)
130
(vector-set! xv i (* 0.9 xv-i)))
133
(define (simplex-shift-list x)
137
(cons (simplex-shift-x x i) (ssl-aux (- i 1)))))
138
(cons x (ssl-aux (- (length x) 1))))
140
; Use the Simplex method to minimize the function (f . x), where
141
; the initial guess is x0 and the fractional tolerance on the value
142
; of the solution is tol.
143
(define (simplex-minimize f x0 tol)
144
(let ((s0 (map (lambda (x) (simplex-point x (apply f x)))
145
(simplex-shift-list x0))))
146
(simplex-iterate f s0 tol)))
added
removed
base/simplex.scm
