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