« back to all changes in this revision
Viewing changes to base/math-utils.scm
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Josselin Mouette
- Date: 2002-04-17 10:36:45 UTC
- Revision ID: james.westby@ubuntu.com-20020417103645-29vomjspk4yf4olw
Tags: upstream-2.1
Import upstream version 2.1
- files added:
- autom4te.cache
- base
- doc
- examples
- utils
added
removed
base/math-utils.scm
1
; libctl: flexible Guile-based control files for scientific software
2
; Copyright (C) 1998, 1999, 2000, 2001, 2002, Steven G. Johnson
3
;
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.
8
;
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.
13
;
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.
18
;
19
; Steven G. Johnson can be contacted at stevenj@alum.mit.edu.
20
21
; ****************************************************************
22
; Miscellaneous math utilities
23
24
; Return the arithmetic sequence (list): start start+step ... (n values)
25
(define (arith-sequence start step n)
26
(define (s x n L) ; tail-recursive helper function
27
(if (= n 0)
28
L
29
(s (binary+ x step) (- n 1) (cons x L))))
30
(reverse (s start n '())))
31
32
; Given a list of numbers, linearly interpolates n values between
33
; each pair of numbers.
34
(define (interpolate n nums)
35
(cons
36
(car nums)
37
(fold-right
38
append '()
39
(map
40
(lambda (x y)
41
(reverse (arith-sequence y (binary/ (binary- x y) (+ n 1)) (+ n 1))))
42
(reverse (cdr (reverse nums))) ; nums w/o last value
43
(cdr nums))))) ; nums w/o first value
44
45
; ****************************************************************
46
; Minimization and root-finding utilities (useful in ctl scripts)
47
48
; The routines are:
49
; minimize: minimize a function of one argument
50
; minimize-multiple: minimize a function of multiple arguments
51
; maximize, maximize-multiple : as above, but maximize
52
; find-root: find the root of a function of one argument
53
; All routines use quadratically convergent methods.
54
55
(define min-arg car)
56
(define min-val cdr)
57
(define max-arg min-arg)
58
(define max-val min-val)
59
60
; One-dimensional minimization (using Brent's method):
61
62
; (minimize f tol) : minimize (f x) with fractional tolerance tol
63
; (minimize f tol guess) : as above, but gives starting guess
64
; (minimize f tol x-min x-max) : as above, but gives range to optimize in
65
; (this is preferred)
66
; All variants return a result that contains both the argument and the
67
; value of the function at its minimum.
68
; (min-arg result) : the argument of the function at its minimum
69
; (min-val result) : the value of the function at its minimum
70
71
(define (minimize f tol . min-max)
72
(define (midpoint a b) (* 0.5 (+ a b)))
73
74
(define (quadratic-min-denom x a b fx fa fb)
75
(magnitude (* 2.0 (- (* (- x a) (- fx fb)) (* (- x b) (- fx fa))))))
76
(define (quadratic-min-num x a b fx fa fb)
77
(let ((den (* 2.0 (- (* (- x a) (- fx fb)) (* (- x b) (- fx fa)))))
78
(num (- (* (- x a) (- x a) (- fx fb))
79
(* (- x b) (- x b) (- fx fa)))))
80
(if (> den 0) (- num) num)))
81
82
(define (tol-scale x) (* tol (+ (magnitude x) 1e-6)))
83
(define (converged? x a b)
84
(<= (magnitude (- x (midpoint a b))) (- (* 2 (tol-scale x)) (* 0.5 (- b a)))))
85
86
(define golden-ratio (* 0.5 (- 3 (sqrt 5))))
87
(define (golden-interpolate x a b)
88
(* golden-ratio (if (>= x (midpoint a b)) (- a x) (- b x))))
89
90
(define (sign x) (if (< x 0) -1 1))
91
92
(define (brent-minimize x a b v w fx fv fw prev-step prev-prev-step)
93
(define (guess-step proposed-step)
94
(let ((step (if (> (magnitude proposed-step) (tol-scale x))
95
proposed-step
96
(* (tol-scale x) (sign proposed-step)))))
97
(let ((u (+ x step)))
98
(let ((fu (f u)))
99
(if (<= fu fx)
100
(if (> u x)
101
(brent-minimize u x b w x fu fw fx step prev-step)
102
(brent-minimize u a x w x fu fw fx step prev-step))
103
(let ((new-a (if (< u x) u a))
104
(new-b (if (< u x) b u)))
105
(if (or (<= fu fw) (= w x))
106
(brent-minimize x new-a new-b w u fx fw fu
107
step prev-step)
108
(if (or (<= fu fv) (= v x) (= v w))
109
(brent-minimize x new-a new-b u w fx fu fw
110
step prev-step)
111
(brent-minimize x new-a new-b v w fx fv fw
112
step prev-step)))))))))
113
114
(if (converged? x a b)
115
(cons x fx)
116
(if (> (magnitude prev-prev-step) (tol-scale x))
117
(let ((p (quadratic-min-num x v w fx fv fw))
118
(q (quadratic-min-denom x v w fx fv fw)))
119
(if (or (>= (magnitude p) (magnitude (* 0.5 q prev-prev-step)))
120
(< p (* q (- a x))) (> p (* q (- b x))))
121
(guess-step (golden-interpolate x a b))
122
(guess-step (/ p q))))
123
(guess-step (golden-interpolate x a b)))))
124
125
(define (bracket-minimum a b c fa fb fc)
126
(if (< fb fc)
127
(list a b c fa fb fc)
128
(let ((u (/ (quadratic-min-num b a c fb fa fc)
129
(max (quadratic-min-denom b a c fb fa fc) 1e-20)))
130
(u-max (+ b (* 100 (- c b)))))
131
(cond
132
((positive? (* (- b u) (- u c)))
133
(let ((fu (f u)))
134
(if (< fu fc)
135
(bracket-minimum b u c fb fu fc)
136
(if (> fu fb)
137
(bracket-minimum a b u fa fb fu)
138
(bracket-minimum b c (+ c (* 1.6 (- c b)))
139
fb fc (f (+ c (* 1.6 (- c b)))))))))
140
((positive? (* (- c u) (- u u-max)))
141
(let ((fu (f u)))
142
(if (< fu fc)
143
(bracket-minimum c u (+ c (* 1.6 (- c b)))
144
fc fu (f (+ c (* 1.6 (- c b)))))
145
(bracket-minimum b c u fb fc fu))))
146
((>= (* (- u u-max) (- u-max c)) 0)
147
(bracket-minimum b c u-max fb fc (f u-max)))
148
(else
149
(bracket-minimum b c (+ c (* 1.6 (- c b)))
150
fb fc (f (+ c (* 1.6 (- c b))))))))))
151
152
(if (= (length min-max) 2)
153
(let ((x-min (first min-max))
154
(x-max (second min-max)))
155
(let ((xm (midpoint x-min x-max)))
156
(let ((fm (f xm)))
157
(brent-minimize xm x-min x-max xm xm fm fm fm 0 0))))
158
(let ((a (if (= (length min-max) 1) (first min-max) 1.0)))
159
(let ((b (if (= a 0) 1.0 0)))
160
(let ((fa (f a)) (fb (f b)))
161
(let ((aa (if (> fb fa) b a))
162
(bb (if (> fb fa) a b))
163
(faa (max fa fb))
164
(fbb (max fa fb)))
165
(let ((bracket
166
(bracket-minimum aa bb (+ bb (* 1.6 (- bb aa)))
167
faa fbb (f (+ bb (* 1.6 (- bb aa)))))))
168
(brent-minimize
169
(second bracket)
170
(min (first bracket) (third bracket))
171
(max (first bracket) (third bracket))
172
(first bracket)
173
(third bracket)
174
(fifth bracket)
175
(fourth bracket)
176
(sixth bracket)
177
0 0))))))))
178
179
; (minimize-multiple f tol arg1 arg2 ... argN) :
180
; Minimize a function f of N arguments, given the fractional tolerance
181
; desired and initial guesses for the arguments.
182
;
183
; (min-arg result) : list of argument values at the minimum
184
; (min-val result) : list of function values at the minimum
185
186
(define (minimize-multiple-expert f tol max-iters fmin guess-args arg-scales)
187
(let ((best-val 1e20) (best-args '()))
188
(subplex
189
(lambda (args)
190
(let ((val (apply f args)))
191
(if (or (null? best-args) (< val best-val))
192
(begin
193
(print "extremization: best so far is "
194
val " at " args "\n")
195
(set! best-val val)
196
(set! best-args args)))
197
val))
198
guess-args tol max-iters
199
(if fmin fmin 0.0) (if fmin true false)
200
arg-scales)))
201
202
(define (minimize-multiple f tol . guess-args)
203
(minimize-multiple-expert f tol 999999999 false guess-args '(0.1)))
204
205
; Yet another alternate multi-dimensional minimization (Simplex algorithm).
206
(define (simplex-minimize-multiple f tol . guess-args)
207
(let ((simplex-result (simplex-minimize f guess-args tol)))
208
(cons (simplex-point-x simplex-result)
209
(simplex-point-val simplex-result))))
210
211
; Alternate multi-dimensional minimization (using Powell's method):
212
; (not the default since it seems to have convergence problems sometimes)
213
214
(define (powell-minimize-multiple f tol . guess-args)
215
(define (create-unit-vector i n)
216
(let ((v (make-vector n 0)))
217
(vector-set! v i 1)
218
v))
219
(define (initial-directions n)
220
(make-initialized-list n (lambda (i) (create-unit-vector i n))))
221
222
(define (v- v1 v2) (vector-map - v1 v2))
223
(define (v+ v1 v2) (vector-map + v1 v2))
224
(define (v* s v) (vector-map (lambda (x) (* s x)) v))
225
(define (v-dot v1 v2) (vector-fold-right + 0 (vector-map * v1 v2)))
226
(define (v-norm v) (sqrt (v-dot v v)))
227
(define (unit-v v) (v* (/ (v-norm v)) v))
228
229
(define (fv v) (apply f (vector->list v)))
230
(define guess-vector (list->vector guess-args))
231
(define (f-dir p0 dir) (lambda (x) (fv (v+ p0 (v* x dir)))))
232
233
(define (minimize-dir p0 dir)
234
(let ((min-result (minimize (f-dir p0 dir) tol)))
235
(cons
236
(v+ p0 (v* (min-arg min-result) dir))
237
(min-val min-result))))
238
239
(define (minimize-dirs p0 dirs)
240
(if (null? dirs)
241
(cons p0 '())
242
(let ((min-result (minimize-dir p0 (car dirs))))
243
(let ((min-results (minimize-dirs (min-arg min-result) (cdr dirs))))
244
(cons (min-arg min-results)
245
(cons (min-val min-result) (min-val min-results)))))))
246
247
(define (replace= val vals els el)
248
(if (null? els) '()
249
(if (= (car vals) val)
250
(cons el (cdr els))
251
(cons (car els) (replace= val (cdr vals) (cdr els) el)))))
252
253
; replace direction where largest decrease occurred:
254
(define (update-dirs decreases dirs p0 p)
255
(replace= (apply max decreases) decreases dirs (v- p p0)))
256
257
(define (minimize-aux p0 fp0 dirs)
258
(let ((min-results (minimize-dirs p0 dirs)))
259
(let ((decreases (map (lambda (val) (- fp0 val)) (min-val min-results)))
260
(p (min-arg min-results))
261
(fp (first (reverse (min-val min-results)))))
262
(if (<= (v-norm (v- p p0))
263
(* tol 0.5 (+ (v-norm p) (v-norm p0) 1e-20)))
264
(cons (vector->list p) fp)
265
(let ((min-result (minimize-dir p (v- p p0))))
266
(minimize-aux (min-arg min-result) (min-val min-result)
267
(update-dirs decreases dirs p0 p)))))))
268
269
(minimize-aux guess-vector (fv guess-vector)
270
(initial-directions (length guess-args))))
271
272
; Maximization variants of the minimize functions:
273
274
(define (maximize f tol . min-max)
275
(let ((result (apply minimize (append (list (compose - f) tol) min-max))))
276
(cons (min-arg result) (- (min-val result)))))
277
278
(define (maximize-multiple f tol . guess-args)
279
(let ((result (apply minimize-multiple
280
(append (list (compose - f) tol) guess-args))))
281
(cons (min-arg result) (- (min-val result)))))
282
283
; Find a root of a function of one argument using Ridder's method.
284
285
; (find-root f tol x-min x-max) : returns the root of the function (f x),
286
; within a fractional tolerance tol. x-min and x-max must bracket the
287
; root; that is, (f x-min) must have a different sign than (f x-max).
288
289
(define (find-root f tol x-min x-max)
290
(define (midpoint a b) (* 0.5 (+ a b)))
291
(define (sign x) (if (< x 0) -1 1))
292
293
(define (best-bracket a b x1 x2 fa fb f1 f2)
294
(if (positive? (* f1 f2))
295
(if (positive? (* fa f1))
296
(list (max x1 x2) b (if (> x1 x2) f1 f2) fb)
297
(list a (min x1 x2) fa (if (< x1 x2) f1 f2)))
298
(if (< x1 x2)
299
(list x1 x2 f1 f2)
300
(list x2 x1 f2 f1))))
301
302
(define (converged? a b x) (< (min (magnitude (- x a)) (magnitude (- x b)))
303
(* tol (magnitude x))))
304
305
; find the root by Ridder's method:
306
(define (ridder a b fa fb)
307
(if (or (= fa 0) (= fb 0))
308
(if (= fa 0) a b)
309
(begin
310
(if (> (* fa fb) 0)
311
(error "x-min and x-max in find-root must bracket the root!"))
312
(let ((m (midpoint a b)))
313
(let ((fm (f m)))
314
(let ((x (+ m (/ (* (- m a) (sign (- fa fb)) fm)
315
(sqrt (- (* fm fm) (* fa fb)))))))
316
(if (or (= fm 0) (converged? a b x))
317
(if (= fm 0) m x)
318
(let ((fx (f x)))
319
(apply ridder (best-bracket a b x m fa fb fx fm))))))))))
320
321
(ridder x-min x-max (f x-min) (f x-max)))
322
323
; ****************************************************************
324
325
; Numerical differentiation:
326
; Compute the numerical derivative of a function f at x, using
327
; Ridder's method of polynomial extrapolation, described e.g. in
328
; Numerical Recipes in C (section 5.7).
329
330
; This is the basic routine, but we wrap it in another interface below
331
; so that dx and tol can be optional arguments.
332
(define (do-derivative f x dx tol)
333
334
; Using Neville's algorithm, compute successively higher-order
335
; extrapolations of the derivative (the "Neville tableau"):
336
(define (deriv-a a0 prev-a fac fac0)
337
(if (null? prev-a)
338
(list a0)
339
(cons a0 (deriv-a (/ (- (* a0 fac) (car prev-a)) (- fac 1))
340
(cdr prev-a) (* fac fac0) fac0))))
341
342
(define (deriv dx df0 err0 prev-a fac0)
343
(let ((a (deriv-a (/ (- (f (+ x dx)) (f (- x dx))) (* 2 dx))
344
prev-a fac0 fac0)))
345
(if (null? prev-a)
346
(deriv (/ dx (sqrt fac0)) (car a) err0 a fac0)
347
(let* ((errs
348
(map max
349
(map magnitude (map - (cdr a) (reverse (cdr (reverse a)))))
350
(map magnitude (map - (cdr a) prev-a))))
351
(errmin (apply min errs))
352
(err (min errmin err0))
353
(df (if (> err err0)
354
df0
355
(cdr (assoc errmin (map cons errs (cdr a)))))))
356
(if (or (< err (* tol df))
357
(> (magnitude (- (car (reverse a)) (car (reverse prev-a))))
358
(* 2 err)))
359
(list df err)
360
(deriv (/ dx (sqrt fac0)) df err a fac0))))))
361
362
(deriv dx 0 1e30 '() 2))
363
364
(define (do-derivative-wrap f x dx-and-tol)
365
(let ((dx (if (> (length dx-and-tol) 0)
366
(car dx-and-tol)
367
(max (magnitude (* x 0.01)) 0.01)))
368
(tol (if (> (length dx-and-tol) 1)
369
(cadr dx-and-tol)
370
0)))
371
(do-derivative f x dx tol)))
372
373
(define derivative-df car)
374
(define derivative-df-err cadr)
375
(define derivative-d2f caddr)
376
(define derivative-d2f-err cadddr)
377
378
(define (derivative f x . dx-and-tol)
379
(do-derivative-wrap f x dx-and-tol))
380
(define (deriv f x . dx-and-tol)
381
(derivative-df (do-derivative-wrap f x dx-and-tol)))
382
383
; Compute both the first and second derivatives at the same time
384
; (using minimal extra function evaluations).
385
(define (derivative2 f x . dx-and-tol)
386
(define f-memo-tab '())
387
(define (f-memo y)
388
(let ((tab-val (assoc y f-memo-tab)))
389
(if tab-val
390
(cdr tab-val)
391
(let ((fy (f y)))
392
(set! f-memo-tab (cons (cons y fy) f-memo-tab))
393
fy))))
394
(define (f-deriv y)
395
(* 2 (/ (- (f-memo y) (f-memo x))
396
(- y x))))
397
(append
398
(do-derivative-wrap f-memo x dx-and-tol)
399
(do-derivative-wrap f-deriv x dx-and-tol)))
400
401
; ****************************************************************
Loggerhead is a web-based interface for Breezy
Version: 2.0.1