40785
|
1 ;; Calculator for GNU Emacs, part II [calc-alg-3.el]
|
|
2 ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
|
|
3 ;; Written by Dave Gillespie, daveg@synaptics.com.
|
|
4
|
|
5 ;; This file is part of GNU Emacs.
|
|
6
|
|
7 ;; GNU Emacs is distributed in the hope that it will be useful,
|
|
8 ;; but WITHOUT ANY WARRANTY. No author or distributor
|
|
9 ;; accepts responsibility to anyone for the consequences of using it
|
|
10 ;; or for whether it serves any particular purpose or works at all,
|
|
11 ;; unless he says so in writing. Refer to the GNU Emacs General Public
|
|
12 ;; License for full details.
|
|
13
|
|
14 ;; Everyone is granted permission to copy, modify and redistribute
|
|
15 ;; GNU Emacs, but only under the conditions described in the
|
|
16 ;; GNU Emacs General Public License. A copy of this license is
|
|
17 ;; supposed to have been given to you along with GNU Emacs so you
|
|
18 ;; can know your rights and responsibilities. It should be in a
|
|
19 ;; file named COPYING. Among other things, the copyright notice
|
|
20 ;; and this notice must be preserved on all copies.
|
|
21
|
|
22
|
|
23
|
|
24 ;; This file is autoloaded from calc-ext.el.
|
|
25 (require 'calc-ext)
|
|
26
|
|
27 (require 'calc-macs)
|
|
28
|
|
29 (defun calc-Need-calc-alg-3 () nil)
|
|
30
|
|
31
|
|
32 (defun calc-find-root (var)
|
|
33 (interactive "sVariable(s) to solve for: ")
|
|
34 (calc-slow-wrapper
|
|
35 (let ((func (if (calc-is-hyperbolic) 'calcFunc-wroot 'calcFunc-root)))
|
|
36 (if (or (equal var "") (equal var "$"))
|
|
37 (calc-enter-result 2 "root" (list func
|
|
38 (calc-top-n 3)
|
|
39 (calc-top-n 1)
|
|
40 (calc-top-n 2)))
|
|
41 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
|
|
42 (not (string-match "\\[" var)))
|
|
43 (math-read-expr (concat "[" var "]"))
|
|
44 (math-read-expr var))))
|
|
45 (if (eq (car-safe var) 'error)
|
|
46 (error "Bad format in expression: %s" (nth 1 var)))
|
|
47 (calc-enter-result 1 "root" (list func
|
|
48 (calc-top-n 2)
|
|
49 var
|
|
50 (calc-top-n 1)))))))
|
|
51 )
|
|
52
|
|
53 (defun calc-find-minimum (var)
|
|
54 (interactive "sVariable(s) to minimize over: ")
|
|
55 (calc-slow-wrapper
|
|
56 (let ((func (if (calc-is-inverse)
|
|
57 (if (calc-is-hyperbolic)
|
|
58 'calcFunc-wmaximize 'calcFunc-maximize)
|
|
59 (if (calc-is-hyperbolic)
|
|
60 'calcFunc-wminimize 'calcFunc-minimize)))
|
|
61 (tag (if (calc-is-inverse) "max" "min")))
|
|
62 (if (or (equal var "") (equal var "$"))
|
|
63 (calc-enter-result 2 tag (list func
|
|
64 (calc-top-n 3)
|
|
65 (calc-top-n 1)
|
|
66 (calc-top-n 2)))
|
|
67 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
|
|
68 (not (string-match "\\[" var)))
|
|
69 (math-read-expr (concat "[" var "]"))
|
|
70 (math-read-expr var))))
|
|
71 (if (eq (car-safe var) 'error)
|
|
72 (error "Bad format in expression: %s" (nth 1 var)))
|
|
73 (calc-enter-result 1 tag (list func
|
|
74 (calc-top-n 2)
|
|
75 var
|
|
76 (calc-top-n 1)))))))
|
|
77 )
|
|
78
|
|
79 (defun calc-find-maximum (var)
|
|
80 (interactive "sVariable to maximize over: ")
|
|
81 (calc-invert-func)
|
|
82 (calc-find-minimum var)
|
|
83 )
|
|
84
|
|
85
|
|
86 (defun calc-poly-interp (arg)
|
|
87 (interactive "P")
|
|
88 (calc-slow-wrapper
|
|
89 (let ((data (calc-top 2)))
|
|
90 (if (or (consp arg) (eq arg 0) (eq arg 2))
|
|
91 (setq data (cons 'vec (calc-top-list 2 2)))
|
|
92 (or (null arg)
|
|
93 (error "Bad prefix argument")))
|
|
94 (if (calc-is-hyperbolic)
|
|
95 (calc-enter-result 1 "rati" (list 'calcFunc-ratint data (calc-top 1)))
|
|
96 (calc-enter-result 1 "poli" (list 'calcFunc-polint data
|
|
97 (calc-top 1))))))
|
|
98 )
|
|
99
|
|
100
|
|
101 (defun calc-curve-fit (arg &optional model coefnames varnames)
|
|
102 (interactive "P")
|
|
103 (calc-slow-wrapper
|
|
104 (setq calc-aborted-prefix nil)
|
|
105 (let ((func (if (calc-is-inverse) 'calcFunc-xfit
|
|
106 (if (calc-is-hyperbolic) 'calcFunc-efit
|
|
107 'calcFunc-fit)))
|
|
108 key (which 0)
|
|
109 n nvars temp data
|
|
110 (homog nil)
|
|
111 (msgs '( "(Press ? for help)"
|
|
112 "1 = linear or multilinear"
|
|
113 "2-9 = polynomial fits; i = interpolating polynomial"
|
|
114 "p = a x^b, ^ = a b^x"
|
|
115 "e = a exp(b x), x = exp(a + b x), l = a + b ln(x)"
|
|
116 "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
|
|
117 "q = a + b (x-c)^2"
|
|
118 "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
|
|
119 "h prefix = homogeneous model (no constant term)"
|
|
120 "' = alg entry, $ = stack, u = Model1, U = Model2")))
|
|
121 (while (not model)
|
|
122 (message "Fit to model: %s:%s"
|
|
123 (nth which msgs)
|
|
124 (if homog " h" ""))
|
|
125 (setq key (read-char))
|
|
126 (cond ((= key ?\C-g)
|
|
127 (keyboard-quit))
|
|
128 ((= key ??)
|
|
129 (setq which (% (1+ which) (length msgs))))
|
|
130 ((memq key '(?h ?H))
|
|
131 (setq homog (not homog)))
|
|
132 ((progn
|
|
133 (if (eq key ?\$)
|
|
134 (setq n 1)
|
|
135 (setq n 0))
|
|
136 (cond ((null arg)
|
|
137 (setq n (1+ n)
|
|
138 data (calc-top n)))
|
|
139 ((or (consp arg) (eq arg 0))
|
|
140 (setq n (+ n 2)
|
|
141 data (calc-top n)
|
|
142 data (if (math-matrixp data)
|
|
143 (append data (list (calc-top (1- n))))
|
|
144 (list 'vec data (calc-top (1- n))))))
|
|
145 ((> (setq arg (prefix-numeric-value arg)) 0)
|
|
146 (setq data (cons 'vec (calc-top-list arg (1+ n)))
|
|
147 n (+ n arg)))
|
|
148 (t (error "Bad prefix argument")))
|
|
149 (or (math-matrixp data) (not (cdr (cdr data)))
|
|
150 (error "Data matrix is not a matrix!"))
|
|
151 (setq nvars (- (length data) 2)
|
|
152 coefnames nil
|
|
153 varnames nil)
|
|
154 nil))
|
|
155 ((= key ?1) ; linear or multilinear
|
|
156 (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
|
|
157 (setq model (math-mul coefnames
|
|
158 (cons 'vec (cons 1 (cdr varnames))))))
|
|
159 ((and (>= key ?2) (<= key ?9)) ; polynomial
|
|
160 (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
|
|
161 (setq model (math-build-polynomial-expr (cdr coefnames)
|
|
162 (nth 1 varnames))))
|
|
163 ((= key ?i) ; exact polynomial
|
|
164 (calc-get-fit-variables 1 (1- (length (nth 1 data)))
|
|
165 (and homog 0))
|
|
166 (setq model (math-build-polynomial-expr (cdr coefnames)
|
|
167 (nth 1 varnames))))
|
|
168 ((= key ?p) ; power law
|
|
169 (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
|
|
170 (setq model (math-mul (nth 1 coefnames)
|
|
171 (calcFunc-reduce
|
|
172 '(var mul var-mul)
|
|
173 (calcFunc-map
|
|
174 '(var pow var-pow)
|
|
175 varnames
|
|
176 (cons 'vec (cdr (cdr coefnames))))))))
|
|
177 ((= key ?^) ; exponential law
|
|
178 (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
|
|
179 (setq model (math-mul (nth 1 coefnames)
|
|
180 (calcFunc-reduce
|
|
181 '(var mul var-mul)
|
|
182 (calcFunc-map
|
|
183 '(var pow var-pow)
|
|
184 (cons 'vec (cdr (cdr coefnames)))
|
|
185 varnames)))))
|
|
186 ((memq key '(?e ?E))
|
|
187 (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
|
|
188 (setq model (math-mul (nth 1 coefnames)
|
|
189 (calcFunc-reduce
|
|
190 '(var mul var-mul)
|
|
191 (calcFunc-map
|
|
192 (if (eq key ?e)
|
|
193 '(var exp var-exp)
|
|
194 '(calcFunc-lambda
|
|
195 (var a var-a)
|
|
196 (^ 10 (var a var-a))))
|
|
197 (calcFunc-map
|
|
198 '(var mul var-mul)
|
|
199 (cons 'vec (cdr (cdr coefnames)))
|
|
200 varnames))))))
|
|
201 ((memq key '(?x ?X))
|
|
202 (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
|
|
203 (setq model (math-mul coefnames
|
|
204 (cons 'vec (cons 1 (cdr varnames)))))
|
|
205 (setq model (if (eq key ?x)
|
|
206 (list 'calcFunc-exp model)
|
|
207 (list '^ 10 model))))
|
|
208 ((memq key '(?l ?L))
|
|
209 (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
|
|
210 (setq model (math-mul coefnames
|
|
211 (cons 'vec
|
|
212 (cons 1 (cdr (calcFunc-map
|
|
213 (if (eq key ?l)
|
|
214 '(var ln var-ln)
|
|
215 '(var log10
|
|
216 var-log10))
|
|
217 varnames)))))))
|
|
218 ((= key ?q)
|
|
219 (calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0))
|
|
220 (let ((c coefnames)
|
|
221 (v varnames))
|
|
222 (setq model (nth 1 c))
|
|
223 (while (setq v (cdr v) c (cdr (cdr c)))
|
|
224 (setq model (math-add
|
|
225 model
|
|
226 (list '*
|
|
227 (car c)
|
|
228 (list '^
|
|
229 (list '- (car v) (nth 1 c))
|
|
230 2)))))))
|
|
231 ((= key ?g)
|
|
232 (setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
|
|
233 varnames '(vec (var XFit var-XFit))
|
|
234 coefnames '(vec (var AFit var-AFit)
|
|
235 (var BFit var-BFit)
|
|
236 (var CFit var-CFit)))
|
|
237 (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))
|
|
238 ((memq key '(?\$ ?\' ?u ?U))
|
|
239 (let* ((defvars nil)
|
|
240 (record-entry nil))
|
|
241 (if (eq key ?\')
|
|
242 (let* ((calc-dollar-values calc-arg-values)
|
|
243 (calc-dollar-used 0)
|
|
244 (calc-hashes-used 0))
|
|
245 (setq model (calc-do-alg-entry "" "Model formula: "))
|
|
246 (if (/= (length model) 1)
|
|
247 (error "Bad format"))
|
|
248 (setq model (car model)
|
|
249 record-entry t)
|
|
250 (if (> calc-dollar-used 0)
|
|
251 (setq coefnames
|
|
252 (cons 'vec
|
|
253 (nthcdr (- (length calc-arg-values)
|
|
254 calc-dollar-used)
|
|
255 (reverse calc-arg-values))))
|
|
256 (if (> calc-hashes-used 0)
|
|
257 (setq coefnames
|
|
258 (cons 'vec (calc-invent-args
|
|
259 calc-hashes-used))))))
|
|
260 (progn
|
|
261 (setq model (cond ((eq key ?u)
|
|
262 (calc-var-value 'var-Model1))
|
|
263 ((eq key ?U)
|
|
264 (calc-var-value 'var-Model2))
|
|
265 (t (calc-top 1))))
|
|
266 (or model (error "User model not yet defined"))
|
|
267 (if (math-vectorp model)
|
|
268 (if (and (memq (length model) '(3 4))
|
|
269 (not (math-objvecp (nth 1 model)))
|
|
270 (math-vectorp (nth 2 model))
|
|
271 (or (null (nth 3 model))
|
|
272 (math-vectorp (nth 3 model))))
|
|
273 (setq varnames (nth 2 model)
|
|
274 coefnames (or (nth 3 model)
|
|
275 (cons 'vec
|
|
276 (math-all-vars-but
|
|
277 model varnames)))
|
|
278 model (nth 1 model))
|
|
279 (error "Incorrect model specifier")))))
|
|
280 (or varnames
|
|
281 (let ((with-y (eq (car-safe model) 'calcFunc-eq)))
|
|
282 (if coefnames
|
|
283 (calc-get-fit-variables (if with-y (1+ nvars) nvars)
|
|
284 (1- (length coefnames))
|
|
285 (math-all-vars-but
|
|
286 model coefnames)
|
|
287 nil with-y)
|
|
288 (let* ((coefs (math-all-vars-but model nil))
|
|
289 (vars nil)
|
|
290 (n (- (length coefs) nvars (if with-y 2 1)))
|
|
291 p)
|
|
292 (if (< n 0)
|
|
293 (error "Not enough variables in model"))
|
|
294 (setq p (nthcdr n coefs))
|
|
295 (setq vars (cdr p))
|
|
296 (setcdr p nil)
|
|
297 (calc-get-fit-variables (if with-y (1+ nvars) nvars)
|
|
298 (length coefs)
|
|
299 vars coefs with-y)))))
|
|
300 (if record-entry
|
|
301 (calc-record (list 'vec model varnames coefnames)
|
|
302 "modl"))))
|
|
303 (t (beep))))
|
|
304 (let ((calc-fit-to-trail t))
|
|
305 (calc-enter-result n (substring (symbol-name func) 9)
|
|
306 (list func model
|
|
307 (if (= (length varnames) 2)
|
|
308 (nth 1 varnames)
|
|
309 varnames)
|
|
310 (if (= (length coefnames) 2)
|
|
311 (nth 1 coefnames)
|
|
312 coefnames)
|
|
313 data))
|
|
314 (if (consp calc-fit-to-trail)
|
|
315 (calc-record (calc-normalize calc-fit-to-trail) "parm")))))
|
|
316 )
|
|
317
|
|
318 (defun calc-invent-independent-variables (n &optional but)
|
|
319 (calc-invent-variables n but '(x y z t) "x")
|
|
320 )
|
|
321
|
|
322 (defun calc-invent-parameter-variables (n &optional but)
|
|
323 (calc-invent-variables n but '(a b c d) "a")
|
|
324 )
|
|
325
|
|
326 (defun calc-invent-variables (num but names base)
|
|
327 (let ((vars nil)
|
|
328 (n num) (nn 0)
|
|
329 var)
|
|
330 (while (and (> n 0) names)
|
|
331 (setq var (math-build-var-name (if (consp names)
|
|
332 (car names)
|
|
333 (concat base (setq nn (1+ nn))))))
|
|
334 (or (math-expr-contains (cons 'vec but) var)
|
|
335 (setq vars (cons var vars)
|
|
336 n (1- n)))
|
|
337 (or (symbolp names) (setq names (cdr names))))
|
|
338 (if (= n 0)
|
|
339 (nreverse vars)
|
|
340 (calc-invent-variables num but t base)))
|
|
341 )
|
|
342
|
|
343 (defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
|
|
344 (or (= nv (if with-y (1+ nvars) nvars))
|
|
345 (error "Wrong number of data vectors for this type of model"))
|
|
346 (if (integerp defv)
|
|
347 (setq homog defv
|
|
348 defv nil))
|
|
349 (if homog
|
|
350 (setq nc (1- nc)))
|
|
351 (or defv
|
|
352 (setq defv (calc-invent-independent-variables nv)))
|
|
353 (or defc
|
|
354 (setq defc (calc-invent-parameter-variables nc defv)))
|
|
355 (let ((vars (read-string (format "Fitting variables: (default %s; %s) "
|
|
356 (mapconcat 'symbol-name
|
|
357 (mapcar (function (lambda (v)
|
|
358 (nth 1 v)))
|
|
359 defv)
|
|
360 ",")
|
|
361 (mapconcat 'symbol-name
|
|
362 (mapcar (function (lambda (v)
|
|
363 (nth 1 v)))
|
|
364 defc)
|
|
365 ","))))
|
|
366 (coefs nil))
|
|
367 (setq vars (if (string-match "\\[" vars)
|
|
368 (math-read-expr vars)
|
|
369 (math-read-expr (concat "[" vars "]"))))
|
|
370 (if (eq (car-safe vars) 'error)
|
|
371 (error "Bad format in expression: %s" (nth 2 vars)))
|
|
372 (or (math-vectorp vars)
|
|
373 (error "Expected a variable or vector of variables"))
|
|
374 (if (equal vars '(vec))
|
|
375 (setq vars (cons 'vec defv)
|
|
376 coefs (cons 'vec defc))
|
|
377 (if (math-vectorp (nth 1 vars))
|
|
378 (if (and (= (length vars) 3)
|
|
379 (math-vectorp (nth 2 vars)))
|
|
380 (setq coefs (nth 2 vars)
|
|
381 vars (nth 1 vars))
|
|
382 (error
|
|
383 "Expected independent variables vector, then parameters vector"))
|
|
384 (setq coefs (cons 'vec defc))))
|
|
385 (or (= nv (1- (length vars)))
|
|
386 (and (not with-y) (= (1+ nv) (1- (length vars))))
|
|
387 (error "Expected %d independent variable%s" nv (if (= nv 1) "" "s")))
|
|
388 (or (= nc (1- (length coefs)))
|
|
389 (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
|
|
390 (if homog
|
|
391 (setq coefs (cons 'vec (cons homog (cdr coefs)))))
|
|
392 (if varnames
|
|
393 (setq model (math-multi-subst model (cdr varnames) (cdr vars))))
|
|
394 (if coefnames
|
|
395 (setq model (math-multi-subst model (cdr coefnames) (cdr coefs))))
|
|
396 (setq varnames vars
|
|
397 coefnames coefs))
|
|
398 )
|
|
399
|
|
400
|
|
401
|
|
402
|
|
403 ;;; The following algorithms are from Numerical Recipes chapter 9.
|
|
404
|
|
405 ;;; "rtnewt" with safety kludges
|
|
406 (defun math-newton-root (expr deriv guess orig-guess limit)
|
|
407 (math-working "newton" guess)
|
|
408 (let* ((var-DUMMY guess)
|
|
409 next dval)
|
|
410 (setq next (math-evaluate-expr expr)
|
|
411 dval (math-evaluate-expr deriv))
|
|
412 (if (and (Math-numberp next)
|
|
413 (Math-numberp dval)
|
|
414 (not (Math-zerop dval)))
|
|
415 (progn
|
|
416 (setq next (math-sub guess (math-div next dval)))
|
|
417 (if (math-nearly-equal guess (setq next (math-float next)))
|
|
418 (progn
|
|
419 (setq var-DUMMY next)
|
|
420 (list 'vec next (math-evaluate-expr expr)))
|
|
421 (if (Math-lessp (math-abs-approx (math-sub next orig-guess))
|
|
422 limit)
|
|
423 (math-newton-root expr deriv next orig-guess limit)
|
|
424 (math-reject-arg next "*Newton's method failed to converge"))))
|
|
425 (math-reject-arg next "*Newton's method encountered a singularity")))
|
|
426 )
|
|
427
|
|
428 ;;; Inspired by "rtsafe"
|
|
429 (defun math-newton-search-root (expr deriv guess vguess ostep oostep
|
|
430 low vlow high vhigh)
|
|
431 (let ((var-DUMMY guess)
|
|
432 (better t)
|
|
433 pos step next vnext)
|
|
434 (if guess
|
|
435 (math-working "newton" (list 'intv 0 low high))
|
|
436 (math-working "bisect" (list 'intv 0 low high))
|
|
437 (setq ostep (math-mul-float (math-sub-float high low)
|
|
438 '(float 5 -1))
|
|
439 guess (math-add-float low ostep)
|
|
440 var-DUMMY guess
|
|
441 vguess (math-evaluate-expr expr))
|
|
442 (or (Math-realp vguess)
|
|
443 (progn
|
|
444 (setq ostep (math-mul-float ostep '(float 6 -1))
|
|
445 guess (math-add-float low ostep)
|
|
446 var-DUMMY guess
|
|
447 vguess (math-evaluate-expr expr))
|
|
448 (or (math-realp vguess)
|
|
449 (progn
|
|
450 (setq ostep (math-mul-float ostep '(float 123456 -5))
|
|
451 guess (math-add-float low ostep)
|
|
452 var-DUMMY guess
|
|
453 vguess nil))))))
|
|
454 (or vguess
|
|
455 (setq vguess (math-evaluate-expr expr)))
|
|
456 (or (Math-realp vguess)
|
|
457 (math-reject-arg guess "*Newton's method encountered a singularity"))
|
|
458 (setq vguess (math-float vguess))
|
|
459 (if (eq (Math-negp vlow) (setq pos (Math-posp vguess)))
|
|
460 (setq high guess
|
|
461 vhigh vguess)
|
|
462 (if (eq (Math-negp vhigh) pos)
|
|
463 (setq low guess
|
|
464 vlow vguess)
|
|
465 (setq better nil)))
|
|
466 (if (or (Math-zerop vguess)
|
|
467 (math-nearly-equal low high))
|
|
468 (list 'vec guess vguess)
|
|
469 (setq step (math-evaluate-expr deriv))
|
|
470 (if (and (Math-realp step)
|
|
471 (not (Math-zerop step))
|
|
472 (setq step (math-div-float vguess (math-float step))
|
|
473 next (math-sub-float guess step))
|
|
474 (not (math-lessp-float high next))
|
|
475 (not (math-lessp-float next low)))
|
|
476 (progn
|
|
477 (setq var-DUMMY next
|
|
478 vnext (math-evaluate-expr expr))
|
|
479 (if (or (Math-zerop vnext)
|
|
480 (math-nearly-equal next guess))
|
|
481 (list 'vec next vnext)
|
|
482 (if (and better
|
|
483 (math-lessp-float (math-abs (or oostep
|
|
484 (math-sub-float
|
|
485 high low)))
|
|
486 (math-abs
|
|
487 (math-mul-float '(float 2 0)
|
|
488 step))))
|
|
489 (math-newton-search-root expr deriv nil nil nil ostep
|
|
490 low vlow high vhigh)
|
|
491 (math-newton-search-root expr deriv next vnext step ostep
|
|
492 low vlow high vhigh))))
|
|
493 (if (or (and (Math-posp vlow) (Math-posp vhigh))
|
|
494 (and (Math-negp vlow) (Math-negp vhigh)))
|
|
495 (math-search-root expr deriv low vlow high vhigh)
|
|
496 (math-newton-search-root expr deriv nil nil nil ostep
|
|
497 low vlow high vhigh)))))
|
|
498 )
|
|
499
|
|
500 ;;; Search for a root in an interval with no overt zero crossing.
|
|
501 (defun math-search-root (expr deriv low vlow high vhigh)
|
|
502 (let (found)
|
|
503 (if root-widen
|
|
504 (let ((iters 0)
|
|
505 (iterlim (if (eq root-widen 'point)
|
|
506 (+ calc-internal-prec 10)
|
|
507 20))
|
|
508 (factor (if (eq root-widen 'point)
|
|
509 '(float 9 0)
|
|
510 '(float 16 -1)))
|
|
511 (prev nil) vprev waslow
|
|
512 diff)
|
|
513 (while (or (and (math-posp vlow) (math-posp vhigh))
|
|
514 (and (math-negp vlow) (math-negp vhigh)))
|
|
515 (math-working "widen" (list 'intv 0 low high))
|
|
516 (if (> (setq iters (1+ iters)) iterlim)
|
|
517 (math-reject-arg (list 'intv 0 low high)
|
|
518 "*Unable to bracket root"))
|
|
519 (if (= iters calc-internal-prec)
|
|
520 (setq factor '(float 16 -1)))
|
|
521 (setq diff (math-mul-float (math-sub-float high low) factor))
|
|
522 (if (Math-zerop diff)
|
|
523 (setq high (calcFunc-incr high 10))
|
|
524 (if (math-lessp-float (math-abs vlow) (math-abs vhigh))
|
|
525 (setq waslow t
|
|
526 prev low
|
|
527 low (math-sub low diff)
|
|
528 var-DUMMY low
|
|
529 vprev vlow
|
|
530 vlow (math-evaluate-expr expr))
|
|
531 (setq waslow nil
|
|
532 prev high
|
|
533 high (math-add high diff)
|
|
534 var-DUMMY high
|
|
535 vprev vhigh
|
|
536 vhigh (math-evaluate-expr expr)))))
|
|
537 (if prev
|
|
538 (if waslow
|
|
539 (setq high prev vhigh vprev)
|
|
540 (setq low prev vlow vprev)))
|
|
541 (setq found t))
|
|
542 (or (Math-realp vlow)
|
|
543 (math-reject-arg vlow 'realp))
|
|
544 (or (Math-realp vhigh)
|
|
545 (math-reject-arg vhigh 'realp))
|
|
546 (let ((xvals (list low high))
|
|
547 (yvals (list vlow vhigh))
|
|
548 (pos (Math-posp vlow))
|
|
549 (levels 0)
|
|
550 (step (math-sub-float high low))
|
|
551 xp yp var-DUMMY)
|
|
552 (while (and (<= (setq levels (1+ levels)) 5)
|
|
553 (not found))
|
|
554 (setq xp xvals
|
|
555 yp yvals
|
|
556 step (math-mul-float step '(float 497 -3)))
|
|
557 (while (and (cdr xp) (not found))
|
|
558 (if (Math-realp (car yp))
|
|
559 (setq low (car xp)
|
|
560 vlow (car yp)))
|
|
561 (setq high (math-add-float (car xp) step)
|
|
562 var-DUMMY high
|
|
563 vhigh (math-evaluate-expr expr))
|
|
564 (math-working "search" high)
|
|
565 (if (and (Math-realp vhigh)
|
|
566 (eq (math-negp vhigh) pos))
|
|
567 (setq found t)
|
|
568 (setcdr xp (cons high (cdr xp)))
|
|
569 (setcdr yp (cons vhigh (cdr yp)))
|
|
570 (setq xp (cdr (cdr xp))
|
|
571 yp (cdr (cdr yp))))))))
|
|
572 (if found
|
|
573 (if (Math-zerop vhigh)
|
|
574 (list 'vec high vhigh)
|
|
575 (if (Math-zerop vlow)
|
|
576 (list 'vec low vlow)
|
|
577 (if deriv
|
|
578 (math-newton-search-root expr deriv nil nil nil nil
|
|
579 low vlow high vhigh)
|
|
580 (math-bisect-root expr low vlow high vhigh))))
|
|
581 (math-reject-arg (list 'intv 3 low high)
|
|
582 "*Unable to find a sign change in this interval")))
|
|
583 )
|
|
584
|
|
585 ;;; "rtbis" (but we should be using Brent's method)
|
|
586 (defun math-bisect-root (expr low vlow high vhigh)
|
|
587 (let ((step (math-sub-float high low))
|
|
588 (pos (Math-posp vhigh))
|
|
589 var-DUMMY
|
|
590 mid vmid)
|
|
591 (while (not (or (math-nearly-equal low
|
|
592 (setq step (math-mul-float
|
|
593 step '(float 5 -1))
|
|
594 mid (math-add-float low step)))
|
|
595 (progn
|
|
596 (setq var-DUMMY mid
|
|
597 vmid (math-evaluate-expr expr))
|
|
598 (Math-zerop vmid))))
|
|
599 (math-working "bisect" mid)
|
|
600 (if (eq (Math-posp vmid) pos)
|
|
601 (setq high mid
|
|
602 vhigh vmid)
|
|
603 (setq low mid
|
|
604 vlow vmid)))
|
|
605 (list 'vec mid vmid))
|
|
606 )
|
|
607
|
|
608 ;;; "mnewt"
|
|
609 (defun math-newton-multi (expr jacob n guess orig-guess limit)
|
|
610 (let ((m -1)
|
|
611 (p guess)
|
|
612 p2 expr-val jacob-val next)
|
|
613 (while (< (setq p (cdr p) m (1+ m)) n)
|
|
614 (set (nth 2 (aref math-root-vars m)) (car p)))
|
|
615 (setq expr-val (math-evaluate-expr expr)
|
|
616 jacob-val (math-evaluate-expr jacob))
|
|
617 (or (and (math-constp expr-val)
|
|
618 (math-constp jacob-val))
|
|
619 (math-reject-arg guess "*Newton's method encountered a singularity"))
|
|
620 (setq next (math-add guess (math-div (math-float (math-neg expr-val))
|
|
621 (math-float jacob-val)))
|
|
622 p guess p2 next)
|
|
623 (math-working "newton" next)
|
|
624 (while (and (setq p (cdr p) p2 (cdr p2))
|
|
625 (math-nearly-equal (car p) (car p2))))
|
|
626 (if p
|
|
627 (if (Math-lessp (math-abs-approx (math-sub next orig-guess))
|
|
628 limit)
|
|
629 (math-newton-multi expr jacob n next orig-guess limit)
|
|
630 (math-reject-arg nil "*Newton's method failed to converge"))
|
|
631 (list 'vec next expr-val)))
|
|
632 )
|
|
633
|
|
634 (defvar math-root-vars [(var DUMMY var-DUMMY)])
|
|
635
|
|
636 (defun math-find-root (expr var guess root-widen)
|
|
637 (if (eq (car-safe expr) 'vec)
|
|
638 (let ((n (1- (length expr)))
|
|
639 (calc-symbolic-mode nil)
|
|
640 (var-DUMMY nil)
|
|
641 (jacob (list 'vec))
|
|
642 p p2 m row)
|
|
643 (or (eq (car-safe var) 'vec)
|
|
644 (math-reject-arg var 'vectorp))
|
|
645 (or (= (length var) (1+ n))
|
|
646 (math-dimension-error))
|
|
647 (setq expr (copy-sequence expr))
|
|
648 (while (>= n (length math-root-vars))
|
|
649 (let ((symb (intern (concat "math-root-v"
|
|
650 (int-to-string
|
|
651 (length math-root-vars))))))
|
|
652 (setq math-root-vars (vconcat math-root-vars
|
|
653 (vector (list 'var symb symb))))))
|
|
654 (setq m -1)
|
|
655 (while (< (setq m (1+ m)) n)
|
|
656 (set (nth 2 (aref math-root-vars m)) nil))
|
|
657 (setq m -1 p var)
|
|
658 (while (setq m (1+ m) p (cdr p))
|
|
659 (or (eq (car-safe (car p)) 'var)
|
|
660 (math-reject-arg var "*Expected a variable"))
|
|
661 (setq p2 expr)
|
|
662 (while (setq p2 (cdr p2))
|
|
663 (setcar p2 (math-expr-subst (car p2) (car p)
|
|
664 (aref math-root-vars m)))))
|
|
665 (or (eq (car-safe guess) 'vec)
|
|
666 (math-reject-arg guess 'vectorp))
|
|
667 (or (= (length guess) (1+ n))
|
|
668 (math-dimension-error))
|
|
669 (setq guess (copy-sequence guess)
|
|
670 p guess)
|
|
671 (while (setq p (cdr p))
|
|
672 (or (Math-numberp (car guess))
|
|
673 (math-reject-arg guess 'numberp))
|
|
674 (setcar p (math-float (car p))))
|
|
675 (setq p expr)
|
|
676 (while (setq p (cdr p))
|
|
677 (if (assq (car-safe (car p)) calc-tweak-eqn-table)
|
|
678 (setcar p (math-sub (nth 1 (car p)) (nth 2 (car p)))))
|
|
679 (setcar p (math-evaluate-expr (car p)))
|
|
680 (setq row (list 'vec)
|
|
681 m -1)
|
|
682 (while (< (setq m (1+ m)) n)
|
|
683 (nconc row (list (math-evaluate-expr
|
|
684 (or (calcFunc-deriv (car p)
|
|
685 (aref math-root-vars m)
|
|
686 nil t)
|
|
687 (math-reject-arg
|
|
688 expr
|
|
689 "*Formulas must be differentiable"))))))
|
|
690 (nconc jacob (list row)))
|
|
691 (setq m (math-abs-approx guess))
|
|
692 (math-newton-multi expr jacob n guess guess
|
|
693 (if (math-zerop m) '(float 1 3) (math-mul m 10))))
|
|
694 (or (eq (car-safe var) 'var)
|
|
695 (math-reject-arg var "*Expected a variable"))
|
|
696 (or (math-expr-contains expr var)
|
|
697 (math-reject-arg expr "*Formula does not contain specified variable"))
|
|
698 (if (assq (car expr) calc-tweak-eqn-table)
|
|
699 (setq expr (math-sub (nth 1 expr) (nth 2 expr))))
|
|
700 (math-with-extra-prec 2
|
|
701 (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
|
|
702 (let* ((calc-symbolic-mode nil)
|
|
703 (var-DUMMY nil)
|
|
704 (expr (math-evaluate-expr expr))
|
|
705 (deriv (calcFunc-deriv expr '(var DUMMY var-DUMMY) nil t))
|
|
706 low high vlow vhigh)
|
|
707 (and deriv (setq deriv (math-evaluate-expr deriv)))
|
|
708 (setq guess (math-float guess))
|
|
709 (if (and (math-numberp guess)
|
|
710 deriv)
|
|
711 (math-newton-root expr deriv guess guess
|
|
712 (if (math-zerop guess) '(float 1 6)
|
|
713 (math-mul (math-abs-approx guess) 100)))
|
|
714 (if (Math-realp guess)
|
|
715 (setq low guess
|
|
716 high guess
|
|
717 var-DUMMY guess
|
|
718 vlow (math-evaluate-expr expr)
|
|
719 vhigh vlow
|
|
720 root-widen 'point)
|
|
721 (if (eq (car guess) 'intv)
|
|
722 (progn
|
|
723 (or (math-constp guess) (math-reject-arg guess 'constp))
|
|
724 (setq low (nth 2 guess)
|
|
725 high (nth 3 guess))
|
|
726 (if (memq (nth 1 guess) '(0 1))
|
|
727 (setq low (calcFunc-incr low 1 high)))
|
|
728 (if (memq (nth 1 guess) '(0 2))
|
|
729 (setq high (calcFunc-incr high -1 low)))
|
|
730 (setq var-DUMMY low
|
|
731 vlow (math-evaluate-expr expr)
|
|
732 var-DUMMY high
|
|
733 vhigh (math-evaluate-expr expr)))
|
|
734 (if (math-complexp guess)
|
|
735 (math-reject-arg "*Complex root finder must have derivative")
|
|
736 (math-reject-arg guess 'realp))))
|
|
737 (if (Math-zerop vlow)
|
|
738 (list 'vec low vlow)
|
|
739 (if (Math-zerop vhigh)
|
|
740 (list 'vec high vhigh)
|
|
741 (if (and deriv (Math-numberp vlow) (Math-numberp vhigh))
|
|
742 (math-newton-search-root expr deriv nil nil nil nil
|
|
743 low vlow high vhigh)
|
|
744 (if (or (and (Math-posp vlow) (Math-posp vhigh))
|
|
745 (and (Math-negp vlow) (Math-negp vhigh))
|
|
746 (not (Math-numberp vlow))
|
|
747 (not (Math-numberp vhigh)))
|
|
748 (math-search-root expr deriv low vlow high vhigh)
|
|
749 (math-bisect-root expr low vlow high vhigh)))))))))
|
|
750 )
|
|
751
|
|
752 (defun calcFunc-root (expr var guess)
|
|
753 (math-find-root expr var guess nil)
|
|
754 )
|
|
755
|
|
756 (defun calcFunc-wroot (expr var guess)
|
|
757 (math-find-root expr var guess t)
|
|
758 )
|
|
759
|
|
760
|
|
761
|
|
762
|
|
763 ;;; The following algorithms come from Numerical Recipes, chapter 10.
|
|
764
|
|
765 (defun math-min-eval (expr a)
|
|
766 (if (Math-vectorp a)
|
|
767 (let ((m -1))
|
|
768 (while (setq m (1+ m) a (cdr a))
|
|
769 (set (nth 2 (aref math-min-vars m)) (car a))))
|
|
770 (setq var-DUMMY a))
|
|
771 (setq a (math-evaluate-expr expr))
|
|
772 (if (Math-ratp a)
|
|
773 (math-float a)
|
|
774 (if (eq (car a) 'float)
|
|
775 a
|
|
776 (math-reject-arg a 'realp)))
|
|
777 )
|
|
778
|
|
779
|
|
780 ;;; A bracket for a minimum is a < b < c where f(b) < f(a) and f(b) < f(c).
|
|
781
|
|
782 ;;; "mnbrak"
|
|
783 (defun math-widen-min (expr a b)
|
|
784 (let ((done nil)
|
|
785 (iters 30)
|
|
786 incr c va vb vc u vu r q ulim bc ba qr)
|
|
787 (or b (setq b (math-mul a '(float 101 -2))))
|
|
788 (setq va (math-min-eval expr a)
|
|
789 vb (math-min-eval expr b))
|
|
790 (if (math-lessp-float va vb)
|
|
791 (setq u a a b b u
|
|
792 vu va va vb vb vu))
|
|
793 (setq c (math-add-float b (math-mul-float '(float 161803 -5)
|
|
794 (math-sub-float b a)))
|
|
795 vc (math-min-eval expr c))
|
|
796 (while (and (not done) (math-lessp-float vc vb))
|
|
797 (math-working "widen" (list 'intv 0 a c))
|
|
798 (if (= (setq iters (1- iters)) 0)
|
|
799 (math-reject-arg nil (format "*Unable to find a %s near the interval"
|
|
800 math-min-or-max)))
|
|
801 (setq bc (math-sub-float b c)
|
|
802 ba (math-sub-float b a)
|
|
803 r (math-mul-float ba (math-sub-float vb vc))
|
|
804 q (math-mul-float bc (math-sub-float vb va))
|
|
805 qr (math-sub-float q r))
|
|
806 (if (math-lessp-float (math-abs qr) '(float 1 -20))
|
|
807 (setq qr (if (math-negp qr) '(float -1 -20) '(float 1 -20))))
|
|
808 (setq u (math-sub-float
|
|
809 b
|
|
810 (math-div-float (math-sub-float (math-mul-float bc q)
|
|
811 (math-mul-float ba r))
|
|
812 (math-mul-float '(float 2 0) qr)))
|
|
813 ulim (math-add-float b (math-mul-float '(float -1 2) bc))
|
|
814 incr (math-negp bc))
|
|
815 (if (if incr (math-lessp-float b u) (math-lessp-float u b))
|
|
816 (if (if incr (math-lessp-float u c) (math-lessp-float c u))
|
|
817 (if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
|
|
818 (setq a b va vb
|
|
819 b u vb vu
|
|
820 done t)
|
|
821 (if (math-lessp-float vb vu)
|
|
822 (setq c u vc vu
|
|
823 done t)
|
|
824 (setq u (math-add-float c (math-mul-float '(float -161803 -5)
|
|
825 bc))
|
|
826 vu (math-min-eval expr u))))
|
|
827 (if (if incr (math-lessp-float u ulim) (math-lessp-float ulim u))
|
|
828 (if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
|
|
829 (setq b c vb vc
|
|
830 c u vc vu
|
|
831 u (math-add-float c (math-mul-float
|
|
832 '(float -161803 -5)
|
|
833 (math-sub-float b c)))
|
|
834 vu (math-min-eval expr u)))
|
|
835 (setq u ulim
|
|
836 vu (math-min-eval expr u))))
|
|
837 (setq u (math-add-float c (math-mul-float '(float -161803 -5)
|
|
838 bc))
|
|
839 vu (math-min-eval expr u)))
|
|
840 (setq a b va vb
|
|
841 b c vb vc
|
|
842 c u vc vu))
|
|
843 (if (math-lessp-float a c)
|
|
844 (list a va b vb c vc)
|
|
845 (list c vc b vb a va)))
|
|
846 )
|
|
847
|
|
848 (defun math-narrow-min (expr a c intv)
|
|
849 (let ((xvals (list a c))
|
|
850 (yvals (list (math-min-eval expr a)
|
|
851 (math-min-eval expr c)))
|
|
852 (levels 0)
|
|
853 (step (math-sub-float c a))
|
|
854 (found nil)
|
|
855 xp yp b)
|
|
856 (while (and (<= (setq levels (1+ levels)) 5)
|
|
857 (not found))
|
|
858 (setq xp xvals
|
|
859 yp yvals
|
|
860 step (math-mul-float step '(float 497 -3)))
|
|
861 (while (and (cdr xp) (not found))
|
|
862 (setq b (math-add-float (car xp) step))
|
|
863 (math-working "search" b)
|
|
864 (setcdr xp (cons b (cdr xp)))
|
|
865 (setcdr yp (cons (math-min-eval expr b) (cdr yp)))
|
|
866 (if (and (math-lessp-float (nth 1 yp) (car yp))
|
|
867 (math-lessp-float (nth 1 yp) (nth 2 yp)))
|
|
868 (setq found t)
|
|
869 (setq xp (cdr xp)
|
|
870 yp (cdr yp))
|
|
871 (if (and (cdr (cdr yp))
|
|
872 (math-lessp-float (nth 1 yp) (car yp))
|
|
873 (math-lessp-float (nth 1 yp) (nth 2 yp)))
|
|
874 (setq found t)
|
|
875 (setq xp (cdr xp)
|
|
876 yp (cdr yp))))))
|
|
877 (if found
|
|
878 (list (car xp) (car yp)
|
|
879 (nth 1 xp) (nth 1 yp)
|
|
880 (nth 2 xp) (nth 2 yp))
|
|
881 (or (if (math-lessp-float (car yvals) (nth 1 yvals))
|
|
882 (and (memq (nth 1 intv) '(2 3))
|
|
883 (let ((min (car yvals)))
|
|
884 (while (and (setq yvals (cdr yvals))
|
|
885 (math-lessp-float min (car yvals))))
|
|
886 (and (not yvals)
|
|
887 (list (nth 2 intv) min))))
|
|
888 (and (memq (nth 1 intv) '(1 3))
|
|
889 (setq yvals (nreverse yvals))
|
|
890 (let ((min (car yvals)))
|
|
891 (while (and (setq yvals (cdr yvals))
|
|
892 (math-lessp-float min (car yvals))))
|
|
893 (and (not yvals)
|
|
894 (list (nth 3 intv) min)))))
|
|
895 (math-reject-arg nil (format "*Unable to find a %s in the interval"
|
|
896 math-min-or-max)))))
|
|
897 )
|
|
898
|
|
899 ;;; "brent"
|
|
900 (defun math-brent-min (expr prec a va x vx b vb)
|
|
901 (let ((iters (+ 20 (* 5 prec)))
|
|
902 (w x)
|
|
903 (vw vx)
|
|
904 (v x)
|
|
905 (vv vx)
|
|
906 (tol (list 'float 1 (- -1 prec)))
|
|
907 (zeps (list 'float 1 (- -5 prec)))
|
|
908 (e '(float 0 0))
|
|
909 u vu xm tol1 tol2 etemp p q r xv xw)
|
|
910 (while (progn
|
|
911 (setq xm (math-mul-float '(float 5 -1)
|
|
912 (math-add-float a b))
|
|
913 tol1 (math-add-float
|
|
914 zeps
|
|
915 (math-mul-float tol (math-abs x)))
|
|
916 tol2 (math-mul-float tol1 '(float 2 0)))
|
|
917 (math-lessp-float (math-sub-float tol2
|
|
918 (math-mul-float
|
|
919 '(float 5 -1)
|
|
920 (math-sub-float b a)))
|
|
921 (math-abs (math-sub-float x xm))))
|
|
922 (if (= (setq iters (1- iters)) 0)
|
|
923 (math-reject-arg nil (format "*Unable to converge on a %s"
|
|
924 math-min-or-max)))
|
|
925 (math-working "brent" x)
|
|
926 (if (math-lessp-float (math-abs e) tol1)
|
|
927 (setq e (if (math-lessp-float x xm)
|
|
928 (math-sub-float b x)
|
|
929 (math-sub-float a x))
|
|
930 d (math-mul-float '(float 381966 -6) e))
|
|
931 (setq xw (math-sub-float x w)
|
|
932 r (math-mul-float xw (math-sub-float vx vv))
|
|
933 xv (math-sub-float x v)
|
|
934 q (math-mul-float xv (math-sub-float vx vw))
|
|
935 p (math-sub-float (math-mul-float xv q)
|
|
936 (math-mul-float xw r))
|
|
937 q (math-mul-float '(float 2 0) (math-sub-float q r)))
|
|
938 (if (math-posp q)
|
|
939 (setq p (math-neg-float p))
|
|
940 (setq q (math-neg-float q)))
|
|
941 (setq etemp e
|
|
942 e d)
|
|
943 (if (and (math-lessp-float (math-abs p)
|
|
944 (math-abs (math-mul-float
|
|
945 '(float 5 -1)
|
|
946 (math-mul-float q etemp))))
|
|
947 (math-lessp-float (math-mul-float
|
|
948 q (math-sub-float a x)) p)
|
|
949 (math-lessp-float p (math-mul-float
|
|
950 q (math-sub-float b x))))
|
|
951 (progn
|
|
952 (setq d (math-div-float p q)
|
|
953 u (math-add-float x d))
|
|
954 (if (or (math-lessp-float (math-sub-float u a) tol2)
|
|
955 (math-lessp-float (math-sub-float b u) tol2))
|
|
956 (setq d (if (math-lessp-float xm x)
|
|
957 (math-neg-float tol1)
|
|
958 tol1))))
|
|
959 (setq e (if (math-lessp-float x xm)
|
|
960 (math-sub-float b x)
|
|
961 (math-sub-float a x))
|
|
962 d (math-mul-float '(float 381966 -6) e))))
|
|
963 (setq u (math-add-float x
|
|
964 (if (math-lessp-float (math-abs d) tol1)
|
|
965 (if (math-negp d)
|
|
966 (math-neg-float tol1)
|
|
967 tol1)
|
|
968 d))
|
|
969 vu (math-min-eval expr u))
|
|
970 (if (math-lessp-float vx vu)
|
|
971 (progn
|
|
972 (if (math-lessp-float u x)
|
|
973 (setq a u)
|
|
974 (setq b u))
|
|
975 (if (or (equal w x)
|
|
976 (not (math-lessp-float vw vu)))
|
|
977 (setq v w vv vw
|
|
978 w u vw vu)
|
|
979 (if (or (equal v x)
|
|
980 (equal v w)
|
|
981 (not (math-lessp-float vv vu)))
|
|
982 (setq v u vv vu))))
|
|
983 (if (math-lessp-float u x)
|
|
984 (setq b x)
|
|
985 (setq a x))
|
|
986 (setq v w vv vw
|
|
987 w x vw vx
|
|
988 x u vx vu)))
|
|
989 (list 'vec x vx))
|
|
990 )
|
|
991
|
|
992 ;;; "powell"
|
|
993 (defun math-powell-min (expr n guesses prec)
|
|
994 (let* ((f1dim (math-line-min-func expr n))
|
|
995 (xi (calcFunc-idn 1 n))
|
|
996 (p (cons 'vec (mapcar 'car guesses)))
|
|
997 (pt p)
|
|
998 (ftol (list 'float 1 (- prec)))
|
|
999 (fret (math-min-eval expr p))
|
|
1000 fp ptt fptt xit i ibig del diff res)
|
|
1001 (while (progn
|
|
1002 (setq fp fret
|
|
1003 ibig 0
|
|
1004 del '(float 0 0)
|
|
1005 i 0)
|
|
1006 (while (<= (setq i (1+ i)) n)
|
|
1007 (setq fptt fret
|
|
1008 res (math-line-min f1dim p
|
|
1009 (math-mat-col xi i)
|
|
1010 n prec)
|
|
1011 p (let ((calc-internal-prec prec))
|
|
1012 (math-normalize (car res)))
|
|
1013 fret (nth 2 res)
|
|
1014 diff (math-abs (math-sub-float fptt fret)))
|
|
1015 (if (math-lessp-float del diff)
|
|
1016 (setq del diff
|
|
1017 ibig i)))
|
|
1018 (math-lessp-float
|
|
1019 (math-mul-float ftol
|
|
1020 (math-add-float (math-abs fp)
|
|
1021 (math-abs fret)))
|
|
1022 (math-mul-float '(float 2 0)
|
|
1023 (math-abs (math-sub-float fp
|
|
1024 fret)))))
|
|
1025 (setq ptt (math-sub (math-mul '(float 2 0) p) pt)
|
|
1026 xit (math-sub p pt)
|
|
1027 pt p
|
|
1028 fptt (math-min-eval expr ptt))
|
|
1029 (if (and (math-lessp-float fptt fp)
|
|
1030 (math-lessp-float
|
|
1031 (math-mul-float
|
|
1032 (math-mul-float '(float 2 0)
|
|
1033 (math-add-float
|
|
1034 (math-sub-float fp
|
|
1035 (math-mul-float '(float 2 0)
|
|
1036 fret))
|
|
1037 fptt))
|
|
1038 (math-sqr-float (math-sub-float
|
|
1039 (math-sub-float fp fret) del)))
|
|
1040 (math-mul-float del
|
|
1041 (math-sqr-float (math-sub-float fp fptt)))))
|
|
1042 (progn
|
|
1043 (setq res (math-line-min f1dim p xit n prec)
|
|
1044 p (car res)
|
|
1045 fret (nth 2 res)
|
|
1046 i 0)
|
|
1047 (while (<= (setq i (1+ i)) n)
|
|
1048 (setcar (nthcdr ibig (nth i xi))
|
|
1049 (nth i (nth 1 res)))))))
|
|
1050 (list 'vec p fret))
|
|
1051 )
|
|
1052
|
|
1053 (defun math-line-min-func (expr n)
|
|
1054 (let ((m -1))
|
|
1055 (while (< (setq m (1+ m)) n)
|
|
1056 (set (nth 2 (aref math-min-vars m))
|
|
1057 (list '+
|
|
1058 (list '*
|
|
1059 '(var DUMMY var-DUMMY)
|
|
1060 (list 'calcFunc-mrow '(var line-xi line-xi) (1+ m)))
|
|
1061 (list 'calcFunc-mrow '(var line-p line-p) (1+ m)))))
|
|
1062 (math-evaluate-expr expr))
|
|
1063 )
|
|
1064
|
|
1065 (defun math-line-min (f1dim line-p line-xi n prec)
|
|
1066 (let* ((var-DUMMY nil)
|
|
1067 (expr (math-evaluate-expr f1dim))
|
|
1068 (params (math-widen-min expr '(float 0 0) '(float 1 0)))
|
|
1069 (res (apply 'math-brent-min expr prec params))
|
|
1070 (xi (math-mul (nth 1 res) line-xi)))
|
|
1071 (list (math-add line-p xi) xi (nth 2 res)))
|
|
1072 )
|
|
1073
|
|
1074
|
|
1075 (defvar math-min-vars [(var DUMMY var-DUMMY)])
|
|
1076
|
|
1077 (defun math-find-minimum (expr var guess min-widen)
|
|
1078 (let* ((calc-symbolic-mode nil)
|
|
1079 (n 0)
|
|
1080 (var-DUMMY nil)
|
|
1081 (isvec (math-vectorp var))
|
|
1082 g guesses)
|
|
1083 (or (math-vectorp var)
|
|
1084 (setq var (list 'vec var)))
|
|
1085 (or (math-vectorp guess)
|
|
1086 (setq guess (list 'vec guess)))
|
|
1087 (or (= (length var) (length guess))
|
|
1088 (math-dimension-error))
|
|
1089 (while (setq var (cdr var) guess (cdr guess))
|
|
1090 (or (eq (car-safe (car var)) 'var)
|
|
1091 (math-reject-arg (car vg) "*Expected a variable"))
|
|
1092 (or (math-expr-contains expr (car var))
|
|
1093 (math-reject-arg (car var)
|
|
1094 "*Formula does not contain specified variable"))
|
|
1095 (while (>= (1+ n) (length math-min-vars))
|
|
1096 (let ((symb (intern (concat "math-min-v"
|
|
1097 (int-to-string
|
|
1098 (length math-min-vars))))))
|
|
1099 (setq math-min-vars (vconcat math-min-vars
|
|
1100 (vector (list 'var symb symb))))))
|
|
1101 (set (nth 2 (aref math-min-vars n)) nil)
|
|
1102 (set (nth 2 (aref math-min-vars (1+ n))) nil)
|
|
1103 (if (math-complexp (car guess))
|
|
1104 (setq expr (math-expr-subst expr
|
|
1105 (car var)
|
|
1106 (list '+ (aref math-min-vars n)
|
|
1107 (list '*
|
|
1108 (aref math-min-vars (1+ n))
|
|
1109 '(cplx 0 1))))
|
|
1110 guesses (let ((g (math-float (math-complex (car guess)))))
|
|
1111 (cons (list (nth 2 g) nil nil)
|
|
1112 (cons (list (nth 1 g) nil nil t)
|
|
1113 guesses)))
|
|
1114 n (+ n 2))
|
|
1115 (setq expr (math-expr-subst expr
|
|
1116 (car var)
|
|
1117 (aref math-min-vars n))
|
|
1118 guesses (cons (if (math-realp (car guess))
|
|
1119 (list (math-float (car guess)) nil nil)
|
|
1120 (if (and (eq (car-safe (car guess)) 'intv)
|
|
1121 (math-constp (car guess)))
|
|
1122 (list (math-mul
|
|
1123 (math-add (nth 2 (car guess))
|
|
1124 (nth 3 (car guess)))
|
|
1125 '(float 5 -1))
|
|
1126 (math-float (nth 2 (car guess)))
|
|
1127 (math-float (nth 3 (car guess)))
|
|
1128 (car guess))
|
|
1129 (math-reject-arg (car guess) 'realp)))
|
|
1130 guesses)
|
|
1131 n (1+ n))))
|
|
1132 (setq guesses (nreverse guesses)
|
|
1133 expr (math-evaluate-expr expr))
|
|
1134 (if (= n 1)
|
|
1135 (let* ((params (if (nth 1 (car guesses))
|
|
1136 (if min-widen
|
|
1137 (math-widen-min expr
|
|
1138 (nth 1 (car guesses))
|
|
1139 (nth 2 (car guesses)))
|
|
1140 (math-narrow-min expr
|
|
1141 (nth 1 (car guesses))
|
|
1142 (nth 2 (car guesses))
|
|
1143 (nth 3 (car guesses))))
|
|
1144 (math-widen-min expr
|
|
1145 (car (car guesses))
|
|
1146 nil)))
|
|
1147 (prec calc-internal-prec)
|
|
1148 (res (if (cdr (cdr params))
|
|
1149 (math-with-extra-prec (+ calc-internal-prec 2)
|
|
1150 (apply 'math-brent-min expr prec params))
|
|
1151 (cons 'vec params))))
|
|
1152 (if isvec
|
|
1153 (list 'vec (list 'vec (nth 1 res)) (nth 2 res))
|
|
1154 res))
|
|
1155 (let* ((prec calc-internal-prec)
|
|
1156 (res (math-with-extra-prec (+ calc-internal-prec 2)
|
|
1157 (math-powell-min expr n guesses prec)))
|
|
1158 (p (nth 1 res))
|
|
1159 (vec (list 'vec)))
|
|
1160 (while (setq p (cdr p))
|
|
1161 (if (nth 3 (car guesses))
|
|
1162 (progn
|
|
1163 (nconc vec (list (math-normalize
|
|
1164 (list 'cplx (car p) (nth 1 p)))))
|
|
1165 (setq p (cdr p)
|
|
1166 guesses (cdr guesses)))
|
|
1167 (nconc vec (list (car p))))
|
|
1168 (setq guesses (cdr guesses)))
|
|
1169 (if isvec
|
|
1170 (list 'vec vec (nth 2 res))
|
|
1171 (list 'vec (nth 1 vec) (nth 2 res))))))
|
|
1172 )
|
|
1173 (setq math-min-or-max "minimum")
|
|
1174
|
|
1175 (defun calcFunc-minimize (expr var guess)
|
|
1176 (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
|
|
1177 (math-min-or-max "minimum"))
|
|
1178 (math-find-minimum (math-normalize expr)
|
|
1179 (math-normalize var)
|
|
1180 (math-normalize guess) nil))
|
|
1181 )
|
|
1182
|
|
1183 (defun calcFunc-wminimize (expr var guess)
|
|
1184 (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
|
|
1185 (math-min-or-max "minimum"))
|
|
1186 (math-find-minimum (math-normalize expr)
|
|
1187 (math-normalize var)
|
|
1188 (math-normalize guess) t))
|
|
1189 )
|
|
1190
|
|
1191 (defun calcFunc-maximize (expr var guess)
|
|
1192 (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
|
|
1193 (math-min-or-max "maximum")
|
|
1194 (res (math-find-minimum (math-normalize (math-neg expr))
|
|
1195 (math-normalize var)
|
|
1196 (math-normalize guess) nil)))
|
|
1197 (list 'vec (nth 1 res) (math-neg (nth 2 res))))
|
|
1198 )
|
|
1199
|
|
1200 (defun calcFunc-wmaximize (expr var guess)
|
|
1201 (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
|
|
1202 (math-min-or-max "maximum")
|
|
1203 (res (math-find-minimum (math-normalize (math-neg expr))
|
|
1204 (math-normalize var)
|
|
1205 (math-normalize guess) t)))
|
|
1206 (list 'vec (nth 1 res) (math-neg (nth 2 res))))
|
|
1207 )
|
|
1208
|
|
1209
|
|
1210
|
|
1211
|
|
1212 ;;; The following algorithms come from Numerical Recipes, chapter 3.
|
|
1213
|
|
1214 (defun calcFunc-polint (data x)
|
|
1215 (or (math-matrixp data) (math-reject-arg data 'matrixp))
|
|
1216 (or (= (length data) 3)
|
|
1217 (math-reject-arg data "*Wrong number of data rows"))
|
|
1218 (or (> (length (nth 1 data)) 2)
|
|
1219 (math-reject-arg data "*Too few data points"))
|
|
1220 (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
|
|
1221 (cons 'vec (mapcar (function (lambda (x) (calcFunc-polint data x)))
|
|
1222 (cdr x)))
|
|
1223 (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
|
|
1224 (math-with-extra-prec 2
|
|
1225 (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
|
|
1226 nil))))
|
|
1227 )
|
|
1228 (put 'calcFunc-polint 'math-expandable t)
|
|
1229
|
|
1230
|
|
1231 (defun calcFunc-ratint (data x)
|
|
1232 (or (math-matrixp data) (math-reject-arg data 'matrixp))
|
|
1233 (or (= (length data) 3)
|
|
1234 (math-reject-arg data "*Wrong number of data rows"))
|
|
1235 (or (> (length (nth 1 data)) 2)
|
|
1236 (math-reject-arg data "*Too few data points"))
|
|
1237 (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
|
|
1238 (cons 'vec (mapcar (function (lambda (x) (calcFunc-ratint data x)))
|
|
1239 (cdr x)))
|
|
1240 (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
|
|
1241 (math-with-extra-prec 2
|
|
1242 (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
|
|
1243 (cdr (cdr (cdr (nth 1 data))))))))
|
|
1244 )
|
|
1245 (put 'calcFunc-ratint 'math-expandable t)
|
|
1246
|
|
1247
|
|
1248 (defun math-poly-interp (xa ya x ratp)
|
|
1249 (let ((n (length xa))
|
|
1250 (dif nil)
|
|
1251 (ns nil)
|
|
1252 (xax nil)
|
|
1253 (c (copy-sequence ya))
|
|
1254 (d (copy-sequence ya))
|
|
1255 (i 0)
|
|
1256 (m 0)
|
|
1257 y dy (xp xa) xpm cp dp temp)
|
|
1258 (while (<= (setq i (1+ i)) n)
|
|
1259 (setq xax (cons (math-sub (car xp) x) xax)
|
|
1260 xp (cdr xp)
|
|
1261 temp (math-abs (car xax)))
|
|
1262 (if (or (null dif) (math-lessp temp dif))
|
|
1263 (setq dif temp
|
|
1264 ns i)))
|
|
1265 (setq xax (nreverse xax)
|
|
1266 ns (1- ns)
|
|
1267 y (nth ns ya))
|
|
1268 (if (math-zerop dif)
|
|
1269 (list y 0)
|
|
1270 (while (< (setq m (1+ m)) n)
|
|
1271 (setq i 0
|
|
1272 xp xax
|
|
1273 xpm (nthcdr m xax)
|
|
1274 cp c
|
|
1275 dp d)
|
|
1276 (while (<= (setq i (1+ i)) (- n m))
|
|
1277 (if ratp
|
|
1278 (let ((t2 (math-div (math-mul (car xp) (car dp)) (car xpm))))
|
|
1279 (setq temp (math-div (math-sub (nth 1 cp) (car dp))
|
|
1280 (math-sub t2 (nth 1 cp))))
|
|
1281 (setcar dp (math-mul (nth 1 cp) temp))
|
|
1282 (setcar cp (math-mul t2 temp)))
|
|
1283 (if (math-equal (car xp) (car xpm))
|
|
1284 (math-reject-arg (cons 'vec xa) "*Duplicate X values"))
|
|
1285 (setq temp (math-div (math-sub (nth 1 cp) (car dp))
|
|
1286 (math-sub (car xp) (car xpm))))
|
|
1287 (setcar dp (math-mul (car xpm) temp))
|
|
1288 (setcar cp (math-mul (car xp) temp)))
|
|
1289 (setq cp (cdr cp)
|
|
1290 dp (cdr dp)
|
|
1291 xp (cdr xp)
|
|
1292 xpm (cdr xpm)))
|
|
1293 (if (< (+ ns ns) (- n m))
|
|
1294 (setq dy (nth ns c))
|
|
1295 (setq ns (1- ns)
|
|
1296 dy (nth ns d)))
|
|
1297 (setq y (math-add y dy)))
|
|
1298 (list y dy)))
|
|
1299 )
|
|
1300
|
|
1301
|
|
1302
|
|
1303 ;;; The following algorithms come from Numerical Recipes, chapter 4.
|
|
1304
|
|
1305 (defun calcFunc-ninteg (expr var lo hi)
|
|
1306 (setq lo (math-evaluate-expr lo)
|
|
1307 hi (math-evaluate-expr hi))
|
|
1308 (or (math-numberp lo) (math-infinitep lo) (math-reject-arg lo 'numberp))
|
|
1309 (or (math-numberp hi) (math-infinitep hi) (math-reject-arg hi 'numberp))
|
|
1310 (if (math-lessp hi lo)
|
|
1311 (math-neg (calcFunc-ninteg expr var hi lo))
|
|
1312 (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
|
|
1313 (let ((var-DUMMY nil)
|
|
1314 (calc-symbolic-mode nil)
|
|
1315 (calc-prefer-frac nil)
|
|
1316 (sum 0))
|
|
1317 (setq expr (math-evaluate-expr expr))
|
|
1318 (if (equal lo '(neg (var inf var-inf)))
|
|
1319 (let ((thi (if (math-lessp hi '(float -2 0))
|
|
1320 hi '(float -2 0))))
|
|
1321 (setq sum (math-ninteg-romberg
|
|
1322 'math-ninteg-midpoint expr
|
|
1323 (math-float lo) (math-float thi) 'inf)
|
|
1324 lo thi)))
|
|
1325 (if (equal hi '(var inf var-inf))
|
|
1326 (let ((tlo (if (math-lessp '(float 2 0) lo)
|
|
1327 lo '(float 2 0))))
|
|
1328 (setq sum (math-add sum
|
|
1329 (math-ninteg-romberg
|
|
1330 'math-ninteg-midpoint expr
|
|
1331 (math-float tlo) (math-float hi) 'inf))
|
|
1332 hi tlo)))
|
|
1333 (or (math-equal lo hi)
|
|
1334 (setq sum (math-add sum
|
|
1335 (math-ninteg-romberg
|
|
1336 'math-ninteg-midpoint expr
|
|
1337 (math-float lo) (math-float hi) nil))))
|
|
1338 sum))
|
|
1339 )
|
|
1340
|
|
1341
|
|
1342 ;;; Open Romberg method; "qromo" in section 4.4.
|
|
1343 (defun math-ninteg-romberg (func expr lo hi mode)
|
|
1344 (let ((curh '(float 1 0))
|
|
1345 (h nil)
|
|
1346 (s nil)
|
|
1347 (j 0)
|
|
1348 (ss nil)
|
|
1349 (prec calc-internal-prec)
|
|
1350 (integ-temp nil))
|
|
1351 (math-with-extra-prec 2
|
|
1352 ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.
|
|
1353 (or (while (and (null ss) (<= (setq j (1+ j)) 8))
|
|
1354 (setq s (nconc s (list (funcall func expr lo hi mode)))
|
|
1355 h (nconc h (list curh)))
|
|
1356 (if (>= j 3)
|
|
1357 (let ((res (math-poly-interp h s '(float 0 0) nil)))
|
|
1358 (if (math-lessp (math-abs (nth 1 res))
|
|
1359 (calcFunc-scf (math-abs (car res))
|
|
1360 (- prec)))
|
|
1361 (setq math-ninteg-convergence j
|
|
1362 ss (car res)))))
|
|
1363 (if (>= j 5)
|
|
1364 (setq s (cdr s)
|
|
1365 h (cdr h)))
|
|
1366 (setq curh (math-div-float curh '(float 9 0))))
|
|
1367 ss
|
|
1368 (math-reject-arg nil (format "*Integral failed to converge")))))
|
|
1369 )
|
|
1370
|
|
1371
|
|
1372 (defun math-ninteg-evaluate (expr x mode)
|
|
1373 (if (eq mode 'inf)
|
|
1374 (setq x (math-div '(float 1 0) x)))
|
|
1375 (let* ((var-DUMMY x)
|
|
1376 (res (math-evaluate-expr expr)))
|
|
1377 (or (Math-numberp res)
|
|
1378 (math-reject-arg res "*Integrand does not evaluate to a number"))
|
|
1379 (if (eq mode 'inf)
|
|
1380 (setq res (math-mul res (math-sqr x))))
|
|
1381 res)
|
|
1382 )
|
|
1383
|
|
1384
|
|
1385 (defun math-ninteg-midpoint (expr lo hi mode) ; uses "integ-temp"
|
|
1386 (if (eq mode 'inf)
|
|
1387 (let ((math-infinite-mode t) temp)
|
|
1388 (setq temp (math-div 1 lo)
|
|
1389 lo (math-div 1 hi)
|
|
1390 hi temp)))
|
|
1391 (if integ-temp
|
|
1392 (let* ((it3 (* 3 (car integ-temp)))
|
|
1393 (math-working-step-2 (* 2 (car integ-temp)))
|
|
1394 (math-working-step 0)
|
|
1395 (range (math-sub hi lo))
|
|
1396 (del (math-div range (math-float it3)))
|
|
1397 (del2 (math-add del del))
|
|
1398 (del3 (math-add del del2))
|
|
1399 (x (math-add lo (math-mul '(float 5 -1) del)))
|
|
1400 (sum '(float 0 0))
|
|
1401 (j 0) temp)
|
|
1402 (while (<= (setq j (1+ j)) (car integ-temp))
|
|
1403 (setq math-working-step (1+ math-working-step)
|
|
1404 temp (math-ninteg-evaluate expr x mode)
|
|
1405 math-working-step (1+ math-working-step)
|
|
1406 sum (math-add sum (math-add temp (math-ninteg-evaluate
|
|
1407 expr (math-add x del2)
|
|
1408 mode)))
|
|
1409 x (math-add x del3)))
|
|
1410 (setq integ-temp (list it3
|
|
1411 (math-add (math-div (nth 1 integ-temp)
|
|
1412 '(float 3 0))
|
|
1413 (math-mul sum del)))))
|
|
1414 (setq integ-temp (list 1 (math-mul
|
|
1415 (math-sub hi lo)
|
|
1416 (math-ninteg-evaluate
|
|
1417 expr
|
|
1418 (math-mul (math-add lo hi) '(float 5 -1))
|
|
1419 mode)))))
|
|
1420 (nth 1 integ-temp)
|
|
1421 )
|
|
1422
|
|
1423
|
|
1424
|
|
1425
|
|
1426
|
|
1427 ;;; The following algorithms come from Numerical Recipes, chapter 14.
|
|
1428
|
|
1429 (setq math-dummy-vars [(var DUMMY var-DUMMY)])
|
|
1430 (setq math-dummy-counter 0)
|
|
1431
|
|
1432 (defun math-dummy-variable ()
|
|
1433 (if (= math-dummy-counter (length math-dummy-vars))
|
|
1434 (let ((symb (intern (format "math-dummy-%d" math-dummy-counter))))
|
|
1435 (setq math-dummy-vars (vconcat math-dummy-vars
|
|
1436 (vector (list 'var symb symb))))))
|
|
1437 (set (nth 2 (aref math-dummy-vars math-dummy-counter)) nil)
|
|
1438 (prog1
|
|
1439 (aref math-dummy-vars math-dummy-counter)
|
|
1440 (setq math-dummy-counter (1+ math-dummy-counter)))
|
|
1441 )
|
|
1442
|
|
1443
|
|
1444
|
|
1445 (defun calcFunc-fit (expr vars &optional coefs data)
|
|
1446 (let ((math-in-fit 10))
|
|
1447 (math-with-extra-prec 2
|
|
1448 (math-general-fit expr vars coefs data nil)))
|
|
1449 )
|
|
1450
|
|
1451 (defun calcFunc-efit (expr vars &optional coefs data)
|
|
1452 (let ((math-in-fit 10))
|
|
1453 (math-with-extra-prec 2
|
|
1454 (math-general-fit expr vars coefs data 'sdev)))
|
|
1455 )
|
|
1456
|
|
1457 (defun calcFunc-xfit (expr vars &optional coefs data)
|
|
1458 (let ((math-in-fit 10))
|
|
1459 (math-with-extra-prec 2
|
|
1460 (math-general-fit expr vars coefs data 'full)))
|
|
1461 )
|
|
1462
|
|
1463 (defun math-general-fit (expr vars coefs data mode)
|
|
1464 (let ((calc-simplify-mode nil)
|
|
1465 (math-dummy-counter math-dummy-counter)
|
|
1466 (math-in-fit 1)
|
|
1467 (extended (eq mode 'full))
|
|
1468 (first-coef math-dummy-counter)
|
|
1469 first-var
|
|
1470 (plain-expr expr)
|
|
1471 orig-expr
|
|
1472 have-sdevs need-chisq chisq
|
|
1473 (x-funcs nil)
|
|
1474 (y-filter nil)
|
|
1475 y-dummy
|
|
1476 (coef-filters nil)
|
|
1477 new-coefs
|
|
1478 (xy-values nil)
|
|
1479 (weights nil)
|
|
1480 (var-YVAL nil) (var-YVALX nil)
|
|
1481 covar beta
|
|
1482 n nn m mm v dummy p)
|
|
1483
|
|
1484 ;; Validate and parse arguments.
|
|
1485 (or data
|
|
1486 (if coefs
|
|
1487 (setq data coefs
|
|
1488 coefs nil)
|
|
1489 (if (math-vectorp expr)
|
|
1490 (if (memq (length expr) '(3 4))
|
|
1491 (setq data vars
|
|
1492 vars (nth 2 expr)
|
|
1493 coefs (nth 3 expr)
|
|
1494 expr (nth 1 expr))
|
|
1495 (math-dimension-error))
|
|
1496 (setq data vars
|
|
1497 vars nil
|
|
1498 coefs nil))))
|
|
1499 (or (math-matrixp data) (math-reject-arg data 'matrixp))
|
|
1500 (setq v (1- (length data))
|
|
1501 n (1- (length (nth 1 data))))
|
|
1502 (or (math-vectorp vars) (null vars)
|
|
1503 (setq vars (list 'vec vars)))
|
|
1504 (or (math-vectorp coefs) (null coefs)
|
|
1505 (setq coefs (list 'vec coefs)))
|
|
1506 (or coefs
|
|
1507 (setq coefs (cons 'vec (math-all-vars-but expr vars))))
|
|
1508 (or vars
|
|
1509 (if (<= (1- (length coefs)) v)
|
|
1510 (math-reject-arg coefs "*Not enough variables in model")
|
|
1511 (setq coefs (copy-sequence coefs))
|
|
1512 (let ((p (nthcdr (- (length coefs) v
|
|
1513 (if (eq (car-safe expr) 'calcFunc-eq) 1 0))
|
|
1514 coefs)))
|
|
1515 (setq vars (cons 'vec (cdr p)))
|
|
1516 (setcdr p nil))))
|
|
1517 (or (= (1- (length vars)) v)
|
|
1518 (= (length vars) v)
|
|
1519 (math-reject-arg vars "*Number of variables does not match data"))
|
|
1520 (setq m (1- (length coefs)))
|
|
1521 (if (< m 1)
|
|
1522 (math-reject-arg coefs "*Need at least one parameter"))
|
|
1523
|
|
1524 ;; Rewrite expr in terms of fitparam and fitvar, make into an equation.
|
|
1525 (setq p coefs)
|
|
1526 (while (setq p (cdr p))
|
|
1527 (or (eq (car-safe (car p)) 'var)
|
|
1528 (math-reject-arg (car p) "*Expected a variable"))
|
|
1529 (setq dummy (math-dummy-variable)
|
|
1530 expr (math-expr-subst expr (car p)
|
|
1531 (list 'calcFunc-fitparam
|
|
1532 (- math-dummy-counter first-coef)))))
|
|
1533 (setq first-var math-dummy-counter
|
|
1534 p vars)
|
|
1535 (while (setq p (cdr p))
|
|
1536 (or (eq (car-safe (car p)) 'var)
|
|
1537 (math-reject-arg (car p) "*Expected a variable"))
|
|
1538 (setq dummy (math-dummy-variable)
|
|
1539 expr (math-expr-subst expr (car p)
|
|
1540 (list 'calcFunc-fitvar
|
|
1541 (- math-dummy-counter first-var)))))
|
|
1542 (if (< math-dummy-counter (+ first-var v))
|
|
1543 (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
|
|
1544 (setq y-dummy dummy
|
|
1545 orig-expr expr)
|
|
1546 (or (eq (car-safe expr) 'calcFunc-eq)
|
|
1547 (setq expr (list 'calcFunc-eq (list 'calcFunc-fitvar v) expr)))
|
|
1548
|
|
1549 (let ((calc-symbolic-mode nil))
|
|
1550
|
|
1551 ;; Apply rewrites to put expr into a linear-like form.
|
|
1552 (setq expr (math-evaluate-expr expr)
|
|
1553 expr (math-rewrite (list 'calcFunc-fitmodel expr)
|
|
1554 '(var FitRules var-FitRules))
|
|
1555 math-in-fit 2
|
|
1556 expr (math-evaluate-expr expr))
|
|
1557 (or (and (eq (car-safe expr) 'calcFunc-fitsystem)
|
|
1558 (= (length expr) 4)
|
|
1559 (math-vectorp (nth 2 expr))
|
|
1560 (math-vectorp (nth 3 expr))
|
|
1561 (> (length (nth 2 expr)) 1)
|
|
1562 (= (length (nth 3 expr)) (1+ m)))
|
|
1563 (math-reject-arg plain-expr "*Model expression is too complex"))
|
|
1564 (setq y-filter (nth 1 expr)
|
|
1565 x-funcs (vconcat (cdr (nth 2 expr)))
|
|
1566 coef-filters (nth 3 expr)
|
|
1567 mm (length x-funcs))
|
|
1568 (if (equal y-filter y-dummy)
|
|
1569 (setq y-filter nil))
|
|
1570
|
|
1571 ;; Build the (square) system of linear equations to be solved.
|
|
1572 (setq beta (cons 'vec (make-list mm 0))
|
|
1573 covar (cons 'vec (mapcar 'copy-sequence (make-list mm beta))))
|
|
1574 (let* ((ptrs (vconcat (cdr data)))
|
|
1575 (isigsq 1)
|
|
1576 (xvals (make-vector mm 0))
|
|
1577 (i 0)
|
|
1578 j k xval yval sigmasqr wt covj covjk covk betaj lud)
|
|
1579 (while (<= (setq i (1+ i)) n)
|
|
1580
|
|
1581 ;; Assign various independent variables for this data point.
|
|
1582 (setq j 0
|
|
1583 sigmasqr nil)
|
|
1584 (while (< j v)
|
|
1585 (aset ptrs j (cdr (aref ptrs j)))
|
|
1586 (setq xval (car (aref ptrs j)))
|
|
1587 (if (= j (1- v))
|
|
1588 (if sigmasqr
|
|
1589 (progn
|
|
1590 (if (eq (car-safe xval) 'sdev)
|
|
1591 (setq sigmasqr (math-add (math-sqr (nth 2 xval))
|
|
1592 sigmasqr)
|
|
1593 xval (nth 1 xval)))
|
|
1594 (if y-filter
|
|
1595 (setq xval (math-make-sdev xval
|
|
1596 (math-sqrt sigmasqr))))))
|
|
1597 (if (eq (car-safe xval) 'sdev)
|
|
1598 (setq sigmasqr (math-add (math-sqr (nth 2 xval))
|
|
1599 (or sigmasqr 0))
|
|
1600 xval (nth 1 xval))))
|
|
1601 (set (nth 2 (aref math-dummy-vars (+ first-var j))) xval)
|
|
1602 (setq j (1+ j)))
|
|
1603
|
|
1604 ;; Compute Y value for this data point.
|
|
1605 (if y-filter
|
|
1606 (setq yval (math-evaluate-expr y-filter))
|
|
1607 (setq yval (symbol-value (nth 2 y-dummy))))
|
|
1608 (if (eq (car-safe yval) 'sdev)
|
|
1609 (setq sigmasqr (math-sqr (nth 2 yval))
|
|
1610 yval (nth 1 yval)))
|
|
1611 (if (= i 1)
|
|
1612 (setq have-sdevs sigmasqr
|
|
1613 need-chisq (or extended
|
|
1614 (and (eq mode 'sdev) (not have-sdevs)))))
|
|
1615 (if have-sdevs
|
|
1616 (if sigmasqr
|
|
1617 (progn
|
|
1618 (setq isigsq (math-div 1 sigmasqr))
|
|
1619 (if need-chisq
|
|
1620 (setq weights (cons isigsq weights))))
|
|
1621 (math-reject-arg yval "*Mixed error forms and plain numbers"))
|
|
1622 (if sigmasqr
|
|
1623 (math-reject-arg yval "*Mixed error forms and plain numbers")))
|
|
1624
|
|
1625 ;; Compute X values for this data point and update covar and beta.
|
|
1626 (if (eq (car-safe xval) 'sdev)
|
|
1627 (set (nth 2 y-dummy) (nth 1 xval)))
|
|
1628 (setq j 0
|
|
1629 covj covar
|
|
1630 betaj beta)
|
|
1631 (while (< j mm)
|
|
1632 (setq wt (math-evaluate-expr (aref x-funcs j)))
|
|
1633 (aset xvals j wt)
|
|
1634 (setq wt (math-mul wt isigsq)
|
|
1635 betaj (cdr betaj)
|
|
1636 covjk (car (setq covj (cdr covj)))
|
|
1637 k 0)
|
|
1638 (while (<= k j)
|
|
1639 (setq covjk (cdr covjk))
|
|
1640 (setcar covjk (math-add (car covjk)
|
|
1641 (math-mul wt (aref xvals k))))
|
|
1642 (setq k (1+ k)))
|
|
1643 (setcar betaj (math-add (car betaj) (math-mul wt yval)))
|
|
1644 (setq j (1+ j)))
|
|
1645 (if need-chisq
|
|
1646 (setq xy-values (cons (append xvals (list yval)) xy-values))))
|
|
1647
|
|
1648 ;; Fill in symmetric half of covar matrix.
|
|
1649 (setq j 0
|
|
1650 covj covar)
|
|
1651 (while (< j (1- mm))
|
|
1652 (setq k j
|
|
1653 j (1+ j)
|
|
1654 covjk (nthcdr j (car (setq covj (cdr covj))))
|
|
1655 covk (nthcdr j covar))
|
|
1656 (while (< (setq k (1+ k)) mm)
|
|
1657 (setq covjk (cdr covjk)
|
|
1658 covk (cdr covk))
|
|
1659 (setcar covjk (nth j (car covk))))))
|
|
1660
|
|
1661 ;; Solve the linear system.
|
|
1662 (if mode
|
|
1663 (progn
|
|
1664 (setq covar (math-matrix-inv-raw covar))
|
|
1665 (if covar
|
|
1666 (setq beta (math-mul covar beta))
|
|
1667 (if (math-zerop (math-abs beta))
|
|
1668 (setq covar (calcFunc-diag 0 (1- (length beta))))
|
|
1669 (math-reject-arg orig-expr "*Singular matrix")))
|
|
1670 (or (math-vectorp covar)
|
|
1671 (setq covar (list 'vec (list 'vec covar)))))
|
|
1672 (setq beta (math-div beta covar)))
|
|
1673
|
|
1674 ;; Compute chi-square statistic if necessary.
|
|
1675 (if need-chisq
|
|
1676 (let (bp xp sum)
|
|
1677 (setq chisq 0)
|
|
1678 (while xy-values
|
|
1679 (setq bp beta
|
|
1680 xp (car xy-values)
|
|
1681 sum 0)
|
|
1682 (while (setq bp (cdr bp))
|
|
1683 (setq sum (math-add sum (math-mul (car bp) (car xp)))
|
|
1684 xp (cdr xp)))
|
|
1685 (setq sum (math-sqr (math-sub (car xp) sum)))
|
|
1686 (if weights (setq sum (math-mul sum (car weights))))
|
|
1687 (setq chisq (math-add chisq sum)
|
|
1688 weights (cdr weights)
|
|
1689 xy-values (cdr xy-values)))))
|
|
1690
|
|
1691 ;; Convert coefficients back into original terms.
|
|
1692 (setq new-coefs (copy-sequence beta))
|
|
1693 (let* ((bp new-coefs)
|
|
1694 (cp covar)
|
|
1695 (sigdat 1)
|
|
1696 (math-in-fit 3)
|
|
1697 (j 0))
|
|
1698 (and mode (not have-sdevs)
|
|
1699 (setq sigdat (if (<= n mm)
|
|
1700 0
|
|
1701 (math-div chisq (- n mm)))))
|
|
1702 (if mode
|
|
1703 (while (setq bp (cdr bp))
|
|
1704 (setcar bp (math-make-sdev
|
|
1705 (car bp)
|
|
1706 (math-sqrt (math-mul (nth (setq j (1+ j))
|
|
1707 (car (setq cp (cdr cp))))
|
|
1708 sigdat))))))
|
|
1709 (setq new-coefs (math-evaluate-expr coef-filters))
|
|
1710 (if calc-fit-to-trail
|
|
1711 (let ((bp new-coefs)
|
|
1712 (cp coefs)
|
|
1713 (vec nil))
|
|
1714 (while (setq bp (cdr bp) cp (cdr cp))
|
|
1715 (setq vec (cons (list 'calcFunc-eq (car cp) (car bp)) vec)))
|
|
1716 (setq calc-fit-to-trail (cons 'vec (nreverse vec)))))))
|
|
1717
|
|
1718 ;; Substitute best-fit coefficients back into original formula.
|
|
1719 (setq expr (math-multi-subst
|
|
1720 orig-expr
|
|
1721 (let ((n v)
|
|
1722 (vec nil))
|
|
1723 (while (>= n 1)
|
|
1724 (setq vec (cons (list 'calcFunc-fitvar n) vec)
|
|
1725 n (1- n)))
|
|
1726 (setq n m)
|
|
1727 (while (>= n 1)
|
|
1728 (setq vec (cons (list 'calcFunc-fitparam n) vec)
|
|
1729 n (1- n)))
|
|
1730 vec)
|
|
1731 (append (cdr new-coefs) (cdr vars))))
|
|
1732
|
|
1733 ;; Package the result.
|
|
1734 (math-normalize
|
|
1735 (if extended
|
|
1736 (list 'vec expr beta covar
|
|
1737 (let ((p coef-filters)
|
|
1738 (n 0))
|
|
1739 (while (and (setq n (1+ n) p (cdr p))
|
|
1740 (eq (car-safe (car p)) 'calcFunc-fitdummy)
|
|
1741 (eq (nth 1 (car p)) n)))
|
|
1742 (if p
|
|
1743 coef-filters
|
|
1744 (list 'vec)))
|
|
1745 chisq
|
|
1746 (if (and have-sdevs (> n mm))
|
|
1747 (list 'calcFunc-utpc chisq (- n mm))
|
|
1748 '(var nan var-nan)))
|
|
1749 expr)))
|
|
1750 )
|
|
1751
|
|
1752 (setq math-in-fit 0)
|
|
1753 (setq calc-fit-to-trail nil)
|
|
1754
|
|
1755 (defun calcFunc-fitvar (x)
|
|
1756 (if (>= math-in-fit 2)
|
|
1757 (progn
|
|
1758 (setq x (aref math-dummy-vars (+ first-var x -1)))
|
|
1759 (or (calc-var-value (nth 2 x)) x))
|
|
1760 (math-reject-arg x))
|
|
1761 )
|
|
1762
|
|
1763 (defun calcFunc-fitparam (x)
|
|
1764 (if (>= math-in-fit 2)
|
|
1765 (progn
|
|
1766 (setq x (aref math-dummy-vars (+ first-coef x -1)))
|
|
1767 (or (calc-var-value (nth 2 x)) x))
|
|
1768 (math-reject-arg x))
|
|
1769 )
|
|
1770
|
|
1771 (defun calcFunc-fitdummy (x)
|
|
1772 (if (= math-in-fit 3)
|
|
1773 (nth x new-coefs)
|
|
1774 (math-reject-arg x))
|
|
1775 )
|
|
1776
|
|
1777 (defun calcFunc-hasfitvars (expr)
|
|
1778 (if (Math-primp expr)
|
|
1779 0
|
|
1780 (if (eq (car expr) 'calcFunc-fitvar)
|
|
1781 (nth 1 expr)
|
|
1782 (apply 'max (mapcar 'calcFunc-hasfitvars (cdr expr)))))
|
|
1783 )
|
|
1784
|
|
1785 (defun calcFunc-hasfitparams (expr)
|
|
1786 (if (Math-primp expr)
|
|
1787 0
|
|
1788 (if (eq (car expr) 'calcFunc-fitparam)
|
|
1789 (nth 1 expr)
|
|
1790 (apply 'max (mapcar 'calcFunc-hasfitparams (cdr expr)))))
|
|
1791 )
|
|
1792
|
|
1793
|
|
1794 (defun math-all-vars-but (expr but)
|
|
1795 (let* ((vars (math-all-vars-in expr))
|
|
1796 (p but))
|
|
1797 (while p
|
|
1798 (setq vars (delq (assoc (car-safe p) vars) vars)
|
|
1799 p (cdr p)))
|
|
1800 (sort (mapcar 'car vars)
|
|
1801 (function (lambda (x y) (string< (nth 1 x) (nth 1 y))))))
|
|
1802 )
|
|
1803
|
|
1804 (defun math-all-vars-in (expr)
|
|
1805 (let ((vars nil)
|
|
1806 found)
|
|
1807 (math-all-vars-rec expr)
|
|
1808 vars)
|
|
1809 )
|
|
1810
|
|
1811 (defun math-all-vars-rec (expr)
|
|
1812 (if (Math-primp expr)
|
|
1813 (if (eq (car-safe expr) 'var)
|
|
1814 (or (math-const-var expr)
|
|
1815 (if (setq found (assoc expr vars))
|
|
1816 (setcdr found (1+ (cdr found)))
|
|
1817 (setq vars (cons (cons expr 1) vars)))))
|
|
1818 (while (setq expr (cdr expr))
|
|
1819 (math-all-vars-rec (car expr))))
|
|
1820 )
|
|
1821
|
|
1822
|
|
1823
|
|
1824
|