3 |
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. |
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. |
4 |
|
|
5 |
;; Author: David Gillespie <daveg@synaptics.com> |
;; Author: David Gillespie <daveg@synaptics.com> |
6 |
;; Maintainers: D. Goel <deego@gnufans.org> |
;; Maintainer: Jay Belanger <belanger@truman.edu> |
|
;; Colin Walters <walters@debian.org> |
|
7 |
|
|
8 |
;; This file is part of GNU Emacs. |
;; This file is part of GNU Emacs. |
9 |
|
|
27 |
;;; Code: |
;;; Code: |
28 |
|
|
29 |
;; This file is autoloaded from calc-ext.el. |
;; This file is autoloaded from calc-ext.el. |
|
(require 'calc-ext) |
|
30 |
|
|
31 |
|
(require 'calc-ext) |
32 |
(require 'calc-macs) |
(require 'calc-macs) |
33 |
|
|
|
(defun calc-Need-calc-alg-3 () nil) |
|
|
|
|
|
|
|
34 |
(defun calc-find-root (var) |
(defun calc-find-root (var) |
35 |
(interactive "sVariable(s) to solve for: ") |
(interactive "sVariable(s) to solve for: ") |
36 |
(calc-slow-wrapper |
(calc-slow-wrapper |
95 |
(calc-enter-result 1 "poli" (list 'calcFunc-polint data |
(calc-enter-result 1 "poli" (list 'calcFunc-polint data |
96 |
(calc-top 1))))))) |
(calc-top 1))))))) |
97 |
|
|
98 |
|
;; The variables calc-curve-nvars, calc-curve-varnames, calc-curve-model and calc-curve-coefnames are local to calc-curve-fit, but are |
99 |
|
;; used by calc-get-fit-variables which is called by calc-curve-fit. |
100 |
|
(defvar calc-curve-nvars) |
101 |
|
(defvar calc-curve-varnames) |
102 |
|
(defvar calc-curve-model) |
103 |
|
(defvar calc-curve-coefnames) |
104 |
|
|
105 |
(defun calc-curve-fit (arg &optional model coefnames varnames) |
(defun calc-curve-fit (arg &optional calc-curve-model |
106 |
|
calc-curve-coefnames calc-curve-varnames) |
107 |
(interactive "P") |
(interactive "P") |
108 |
(calc-slow-wrapper |
(calc-slow-wrapper |
109 |
(setq calc-aborted-prefix nil) |
(setq calc-aborted-prefix nil) |
111 |
(if (calc-is-hyperbolic) 'calcFunc-efit |
(if (calc-is-hyperbolic) 'calcFunc-efit |
112 |
'calcFunc-fit))) |
'calcFunc-fit))) |
113 |
key (which 0) |
key (which 0) |
114 |
n nvars temp data |
n calc-curve-nvars temp data |
115 |
(homog nil) |
(homog nil) |
116 |
(msgs '( "(Press ? for help)" |
(msgs '( "(Press ? for help)" |
117 |
"1 = linear or multilinear" |
"1 = linear or multilinear" |
123 |
"g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)" |
"g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)" |
124 |
"h prefix = homogeneous model (no constant term)" |
"h prefix = homogeneous model (no constant term)" |
125 |
"' = alg entry, $ = stack, u = Model1, U = Model2"))) |
"' = alg entry, $ = stack, u = Model1, U = Model2"))) |
126 |
(while (not model) |
(while (not calc-curve-model) |
127 |
(message "Fit to model: %s:%s" |
(message "Fit to model: %s:%s" |
128 |
(nth which msgs) |
(nth which msgs) |
129 |
(if homog " h" "")) |
(if homog " h" "")) |
153 |
(t (error "Bad prefix argument"))) |
(t (error "Bad prefix argument"))) |
154 |
(or (math-matrixp data) (not (cdr (cdr data))) |
(or (math-matrixp data) (not (cdr (cdr data))) |
155 |
(error "Data matrix is not a matrix!")) |
(error "Data matrix is not a matrix!")) |
156 |
(setq nvars (- (length data) 2) |
(setq calc-curve-nvars (- (length data) 2) |
157 |
coefnames nil |
calc-curve-coefnames nil |
158 |
varnames nil) |
calc-curve-varnames nil) |
159 |
nil)) |
nil)) |
160 |
((= key ?1) ; linear or multilinear |
((= key ?1) ; linear or multilinear |
161 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 0)) |
(calc-get-fit-variables calc-curve-nvars |
162 |
(setq model (math-mul coefnames |
(1+ calc-curve-nvars) (and homog 0)) |
163 |
(cons 'vec (cons 1 (cdr varnames)))))) |
(setq calc-curve-model (math-mul calc-curve-coefnames |
164 |
|
(cons 'vec (cons 1 (cdr calc-curve-varnames)))))) |
165 |
((and (>= key ?2) (<= key ?9)) ; polynomial |
((and (>= key ?2) (<= key ?9)) ; polynomial |
166 |
(calc-get-fit-variables 1 (- key ?0 -1) (and homog 0)) |
(calc-get-fit-variables 1 (- key ?0 -1) (and homog 0)) |
167 |
(setq model (math-build-polynomial-expr (cdr coefnames) |
(setq calc-curve-model |
168 |
(nth 1 varnames)))) |
(math-build-polynomial-expr (cdr calc-curve-coefnames) |
169 |
|
(nth 1 calc-curve-varnames)))) |
170 |
((= key ?i) ; exact polynomial |
((= key ?i) ; exact polynomial |
171 |
(calc-get-fit-variables 1 (1- (length (nth 1 data))) |
(calc-get-fit-variables 1 (1- (length (nth 1 data))) |
172 |
(and homog 0)) |
(and homog 0)) |
173 |
(setq model (math-build-polynomial-expr (cdr coefnames) |
(setq calc-curve-model |
174 |
(nth 1 varnames)))) |
(math-build-polynomial-expr (cdr calc-curve-coefnames) |
175 |
|
(nth 1 calc-curve-varnames)))) |
176 |
((= key ?p) ; power law |
((= key ?p) ; power law |
177 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 1)) |
(calc-get-fit-variables calc-curve-nvars |
178 |
(setq model (math-mul (nth 1 coefnames) |
(1+ calc-curve-nvars) (and homog 1)) |
179 |
|
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) |
180 |
(calcFunc-reduce |
(calcFunc-reduce |
181 |
'(var mul var-mul) |
'(var mul var-mul) |
182 |
(calcFunc-map |
(calcFunc-map |
183 |
'(var pow var-pow) |
'(var pow var-pow) |
184 |
varnames |
calc-curve-varnames |
185 |
(cons 'vec (cdr (cdr coefnames)))))))) |
(cons 'vec (cdr (cdr calc-curve-coefnames)))))))) |
186 |
((= key ?^) ; exponential law |
((= key ?^) ; exponential law |
187 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 1)) |
(calc-get-fit-variables calc-curve-nvars |
188 |
(setq model (math-mul (nth 1 coefnames) |
(1+ calc-curve-nvars) (and homog 1)) |
189 |
|
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) |
190 |
(calcFunc-reduce |
(calcFunc-reduce |
191 |
'(var mul var-mul) |
'(var mul var-mul) |
192 |
(calcFunc-map |
(calcFunc-map |
193 |
'(var pow var-pow) |
'(var pow var-pow) |
194 |
(cons 'vec (cdr (cdr coefnames))) |
(cons 'vec (cdr (cdr calc-curve-coefnames))) |
195 |
varnames))))) |
calc-curve-varnames))))) |
196 |
((memq key '(?e ?E)) |
((memq key '(?e ?E)) |
197 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 1)) |
(calc-get-fit-variables calc-curve-nvars |
198 |
(setq model (math-mul (nth 1 coefnames) |
(1+ calc-curve-nvars) (and homog 1)) |
199 |
|
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames) |
200 |
(calcFunc-reduce |
(calcFunc-reduce |
201 |
'(var mul var-mul) |
'(var mul var-mul) |
202 |
(calcFunc-map |
(calcFunc-map |
207 |
(^ 10 (var a var-a)))) |
(^ 10 (var a var-a)))) |
208 |
(calcFunc-map |
(calcFunc-map |
209 |
'(var mul var-mul) |
'(var mul var-mul) |
210 |
(cons 'vec (cdr (cdr coefnames))) |
(cons 'vec (cdr (cdr calc-curve-coefnames))) |
211 |
varnames)))))) |
calc-curve-varnames)))))) |
212 |
((memq key '(?x ?X)) |
((memq key '(?x ?X)) |
213 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 0)) |
(calc-get-fit-variables calc-curve-nvars |
214 |
(setq model (math-mul coefnames |
(1+ calc-curve-nvars) (and homog 0)) |
215 |
(cons 'vec (cons 1 (cdr varnames))))) |
(setq calc-curve-model (math-mul calc-curve-coefnames |
216 |
(setq model (if (eq key ?x) |
(cons 'vec (cons 1 (cdr calc-curve-varnames))))) |
217 |
(list 'calcFunc-exp model) |
(setq calc-curve-model (if (eq key ?x) |
218 |
(list '^ 10 model)))) |
(list 'calcFunc-exp calc-curve-model) |
219 |
|
(list '^ 10 calc-curve-model)))) |
220 |
((memq key '(?l ?L)) |
((memq key '(?l ?L)) |
221 |
(calc-get-fit-variables nvars (1+ nvars) (and homog 0)) |
(calc-get-fit-variables calc-curve-nvars |
222 |
(setq model (math-mul coefnames |
(1+ calc-curve-nvars) (and homog 0)) |
223 |
|
(setq calc-curve-model (math-mul calc-curve-coefnames |
224 |
(cons 'vec |
(cons 'vec |
225 |
(cons 1 (cdr (calcFunc-map |
(cons 1 (cdr (calcFunc-map |
226 |
(if (eq key ?l) |
(if (eq key ?l) |
227 |
'(var ln var-ln) |
'(var ln var-ln) |
228 |
'(var log10 |
'(var log10 |
229 |
var-log10)) |
var-log10)) |
230 |
varnames))))))) |
calc-curve-varnames))))))) |
231 |
((= key ?q) |
((= key ?q) |
232 |
(calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0)) |
(calc-get-fit-variables calc-curve-nvars |
233 |
(let ((c coefnames) |
(1+ (* 2 calc-curve-nvars)) (and homog 0)) |
234 |
(v varnames)) |
(let ((c calc-curve-coefnames) |
235 |
(setq model (nth 1 c)) |
(v calc-curve-varnames)) |
236 |
|
(setq calc-curve-model (nth 1 c)) |
237 |
(while (setq v (cdr v) c (cdr (cdr c))) |
(while (setq v (cdr v) c (cdr (cdr c))) |
238 |
(setq model (math-add |
(setq calc-curve-model (math-add |
239 |
model |
calc-curve-model |
240 |
(list '* |
(list '* |
241 |
(car c) |
(car c) |
242 |
(list '^ |
(list '^ |
243 |
(list '- (car v) (nth 1 c)) |
(list '- (car v) (nth 1 c)) |
244 |
2))))))) |
2))))))) |
245 |
((= key ?g) |
((= key ?g) |
246 |
(setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") |
(setq calc-curve-model |
247 |
varnames '(vec (var XFit var-XFit)) |
(math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") |
248 |
coefnames '(vec (var AFit var-AFit) |
calc-curve-varnames '(vec (var XFit var-XFit)) |
249 |
|
calc-curve-coefnames '(vec (var AFit var-AFit) |
250 |
(var BFit var-BFit) |
(var BFit var-BFit) |
251 |
(var CFit var-CFit))) |
(var CFit var-CFit))) |
252 |
(calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1))) |
(calc-get-fit-variables 1 (1- (length calc-curve-coefnames)) |
253 |
|
(and homog 1))) |
254 |
((memq key '(?\$ ?\' ?u ?U)) |
((memq key '(?\$ ?\' ?u ?U)) |
255 |
(let* ((defvars nil) |
(let* ((defvars nil) |
256 |
(record-entry nil)) |
(record-entry nil)) |
258 |
(let* ((calc-dollar-values calc-arg-values) |
(let* ((calc-dollar-values calc-arg-values) |
259 |
(calc-dollar-used 0) |
(calc-dollar-used 0) |
260 |
(calc-hashes-used 0)) |
(calc-hashes-used 0)) |
261 |
(setq model (calc-do-alg-entry "" "Model formula: ")) |
(setq calc-curve-model (calc-do-alg-entry "" "Model formula: ")) |
262 |
(if (/= (length model) 1) |
(if (/= (length calc-curve-model) 1) |
263 |
(error "Bad format")) |
(error "Bad format")) |
264 |
(setq model (car model) |
(setq calc-curve-model (car calc-curve-model) |
265 |
record-entry t) |
record-entry t) |
266 |
(if (> calc-dollar-used 0) |
(if (> calc-dollar-used 0) |
267 |
(setq coefnames |
(setq calc-curve-coefnames |
268 |
(cons 'vec |
(cons 'vec |
269 |
(nthcdr (- (length calc-arg-values) |
(nthcdr (- (length calc-arg-values) |
270 |
calc-dollar-used) |
calc-dollar-used) |
271 |
(reverse calc-arg-values)))) |
(reverse calc-arg-values)))) |
272 |
(if (> calc-hashes-used 0) |
(if (> calc-hashes-used 0) |
273 |
(setq coefnames |
(setq calc-curve-coefnames |
274 |
(cons 'vec (calc-invent-args |
(cons 'vec (calc-invent-args |
275 |
calc-hashes-used)))))) |
calc-hashes-used)))))) |
276 |
(progn |
(progn |
277 |
(setq model (cond ((eq key ?u) |
(setq calc-curve-model (cond ((eq key ?u) |
278 |
(calc-var-value 'var-Model1)) |
(calc-var-value 'var-Model1)) |
279 |
((eq key ?U) |
((eq key ?U) |
280 |
(calc-var-value 'var-Model2)) |
(calc-var-value 'var-Model2)) |
281 |
(t (calc-top 1)))) |
(t (calc-top 1)))) |
282 |
(or model (error "User model not yet defined")) |
(or calc-curve-model (error "User model not yet defined")) |
283 |
(if (math-vectorp model) |
(if (math-vectorp calc-curve-model) |
284 |
(if (and (memq (length model) '(3 4)) |
(if (and (memq (length calc-curve-model) '(3 4)) |
285 |
(not (math-objvecp (nth 1 model))) |
(not (math-objvecp (nth 1 calc-curve-model))) |
286 |
(math-vectorp (nth 2 model)) |
(math-vectorp (nth 2 calc-curve-model)) |
287 |
(or (null (nth 3 model)) |
(or (null (nth 3 calc-curve-model)) |
288 |
(math-vectorp (nth 3 model)))) |
(math-vectorp (nth 3 calc-curve-model)))) |
289 |
(setq varnames (nth 2 model) |
(setq calc-curve-varnames (nth 2 calc-curve-model) |
290 |
coefnames (or (nth 3 model) |
calc-curve-coefnames |
291 |
(cons 'vec |
(or (nth 3 calc-curve-model) |
292 |
(math-all-vars-but |
(cons 'vec |
293 |
model varnames))) |
(math-all-vars-but |
294 |
model (nth 1 model)) |
calc-curve-model calc-curve-varnames))) |
295 |
|
calc-curve-model (nth 1 calc-curve-model)) |
296 |
(error "Incorrect model specifier"))))) |
(error "Incorrect model specifier"))))) |
297 |
(or varnames |
(or calc-curve-varnames |
298 |
(let ((with-y (eq (car-safe model) 'calcFunc-eq))) |
(let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq))) |
299 |
(if coefnames |
(if calc-curve-coefnames |
300 |
(calc-get-fit-variables (if with-y (1+ nvars) nvars) |
(calc-get-fit-variables |
301 |
(1- (length coefnames)) |
(if with-y (1+ calc-curve-nvars) calc-curve-nvars) |
302 |
(math-all-vars-but |
(1- (length calc-curve-coefnames)) |
303 |
model coefnames) |
(math-all-vars-but |
304 |
nil with-y) |
calc-curve-model calc-curve-coefnames) |
305 |
(let* ((coefs (math-all-vars-but model nil)) |
nil with-y) |
306 |
|
(let* ((coefs (math-all-vars-but calc-curve-model nil)) |
307 |
(vars nil) |
(vars nil) |
308 |
(n (- (length coefs) nvars (if with-y 2 1))) |
(n (- (length coefs) calc-curve-nvars (if with-y 2 1))) |
309 |
p) |
p) |
310 |
(if (< n 0) |
(if (< n 0) |
311 |
(error "Not enough variables in model")) |
(error "Not enough variables in model")) |
312 |
(setq p (nthcdr n coefs)) |
(setq p (nthcdr n coefs)) |
313 |
(setq vars (cdr p)) |
(setq vars (cdr p)) |
314 |
(setcdr p nil) |
(setcdr p nil) |
315 |
(calc-get-fit-variables (if with-y (1+ nvars) nvars) |
(calc-get-fit-variables |
316 |
(length coefs) |
(if with-y (1+ calc-curve-nvars) calc-curve-nvars) |
317 |
vars coefs with-y))))) |
(length coefs) |
318 |
|
vars coefs with-y))))) |
319 |
(if record-entry |
(if record-entry |
320 |
(calc-record (list 'vec model varnames coefnames) |
(calc-record (list 'vec calc-curve-model |
321 |
|
calc-curve-varnames calc-curve-coefnames) |
322 |
"modl")))) |
"modl")))) |
323 |
(t (beep)))) |
(t (beep)))) |
324 |
(let ((calc-fit-to-trail t)) |
(let ((calc-fit-to-trail t)) |
325 |
(calc-enter-result n (substring (symbol-name func) 9) |
(calc-enter-result n (substring (symbol-name func) 9) |
326 |
(list func model |
(list func calc-curve-model |
327 |
(if (= (length varnames) 2) |
(if (= (length calc-curve-varnames) 2) |
328 |
(nth 1 varnames) |
(nth 1 calc-curve-varnames) |
329 |
varnames) |
calc-curve-varnames) |
330 |
(if (= (length coefnames) 2) |
(if (= (length calc-curve-coefnames) 2) |
331 |
(nth 1 coefnames) |
(nth 1 calc-curve-coefnames) |
332 |
coefnames) |
calc-curve-coefnames) |
333 |
data)) |
data)) |
334 |
(if (consp calc-fit-to-trail) |
(if (consp calc-fit-to-trail) |
335 |
(calc-record (calc-normalize calc-fit-to-trail) "parm")))))) |
(calc-record (calc-normalize calc-fit-to-trail) "parm")))))) |
358 |
(calc-invent-variables num but t base)))) |
(calc-invent-variables num but t base)))) |
359 |
|
|
360 |
(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog) |
(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog) |
361 |
(or (= nv (if with-y (1+ nvars) nvars)) |
(or (= nv (if with-y (1+ calc-curve-nvars) calc-curve-nvars)) |
362 |
(error "Wrong number of data vectors for this type of model")) |
(error "Wrong number of data vectors for this type of model")) |
363 |
(if (integerp defv) |
(if (integerp defv) |
364 |
(setq homog defv |
(setq homog defv |
406 |
(error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s"))) |
(error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s"))) |
407 |
(if homog |
(if homog |
408 |
(setq coefs (cons 'vec (cons homog (cdr coefs))))) |
(setq coefs (cons 'vec (cons homog (cdr coefs))))) |
409 |
(if varnames |
(if calc-curve-varnames |
410 |
(setq model (math-multi-subst model (cdr varnames) (cdr vars)))) |
(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-varnames) (cdr vars)))) |
411 |
(if coefnames |
(if calc-curve-coefnames |
412 |
(setq model (math-multi-subst model (cdr coefnames) (cdr coefs)))) |
(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-coefnames) (cdr coefs)))) |
413 |
(setq varnames vars |
(setq calc-curve-varnames vars |
414 |
coefnames coefs))) |
calc-curve-coefnames coefs))) |
415 |
|
|
416 |
|
|
417 |
|
|
419 |
;;; The following algorithms are from Numerical Recipes chapter 9. |
;;; The following algorithms are from Numerical Recipes chapter 9. |
420 |
|
|
421 |
;;; "rtnewt" with safety kludges |
;;; "rtnewt" with safety kludges |
422 |
|
|
423 |
|
(defvar var-DUMMY) |
424 |
|
|
425 |
(defun math-newton-root (expr deriv guess orig-guess limit) |
(defun math-newton-root (expr deriv guess orig-guess limit) |
426 |
(math-working "newton" guess) |
(math-working "newton" guess) |
427 |
(let* ((var-DUMMY guess) |
(let* ((var-DUMMY guess) |
515 |
low vlow high vhigh)))))) |
low vlow high vhigh)))))) |
516 |
|
|
517 |
;;; Search for a root in an interval with no overt zero crossing. |
;;; Search for a root in an interval with no overt zero crossing. |
518 |
|
|
519 |
|
;; The variable math-root-widen is local to math-find-root, but |
520 |
|
;; is used by math-search-root, which is called (directly and |
521 |
|
;; indirectly) by math-find-root. |
522 |
|
(defvar math-root-widen) |
523 |
|
|
524 |
(defun math-search-root (expr deriv low vlow high vhigh) |
(defun math-search-root (expr deriv low vlow high vhigh) |
525 |
(let (found) |
(let (found) |
526 |
(if root-widen |
(if math-root-widen |
527 |
(let ((iters 0) |
(let ((iters 0) |
528 |
(iterlim (if (eq root-widen 'point) |
(iterlim (if (eq math-root-widen 'point) |
529 |
(+ calc-internal-prec 10) |
(+ calc-internal-prec 10) |
530 |
20)) |
20)) |
531 |
(factor (if (eq root-widen 'point) |
(factor (if (eq math-root-widen 'point) |
532 |
'(float 9 0) |
'(float 9 0) |
533 |
'(float 16 -1))) |
'(float 16 -1))) |
534 |
(prev nil) vprev waslow |
(prev nil) vprev waslow |
627 |
(list 'vec mid vmid))) |
(list 'vec mid vmid))) |
628 |
|
|
629 |
;;; "mnewt" |
;;; "mnewt" |
630 |
|
|
631 |
|
(defvar math-root-vars [(var DUMMY var-DUMMY)]) |
632 |
|
|
633 |
(defun math-newton-multi (expr jacob n guess orig-guess limit) |
(defun math-newton-multi (expr jacob n guess orig-guess limit) |
634 |
(let ((m -1) |
(let ((m -1) |
635 |
(p guess) |
(p guess) |
654 |
(math-reject-arg nil "*Newton's method failed to converge")) |
(math-reject-arg nil "*Newton's method failed to converge")) |
655 |
(list 'vec next expr-val)))) |
(list 'vec next expr-val)))) |
656 |
|
|
|
(defvar math-root-vars [(var DUMMY var-DUMMY)]) |
|
657 |
|
|
658 |
(defun math-find-root (expr var guess root-widen) |
(defun math-find-root (expr var guess math-root-widen) |
659 |
(if (eq (car-safe expr) 'vec) |
(if (eq (car-safe expr) 'vec) |
660 |
(let ((n (1- (length expr))) |
(let ((n (1- (length expr))) |
661 |
(calc-symbolic-mode nil) |
(calc-symbolic-mode nil) |
739 |
var-DUMMY guess |
var-DUMMY guess |
740 |
vlow (math-evaluate-expr expr) |
vlow (math-evaluate-expr expr) |
741 |
vhigh vlow |
vhigh vlow |
742 |
root-widen 'point) |
math-root-widen 'point) |
743 |
(if (eq (car guess) 'intv) |
(if (eq (car guess) 'intv) |
744 |
(progn |
(progn |
745 |
(or (math-constp guess) (math-reject-arg guess 'constp)) |
(or (math-constp guess) (math-reject-arg guess 'constp)) |
781 |
|
|
782 |
;;; The following algorithms come from Numerical Recipes, chapter 10. |
;;; The following algorithms come from Numerical Recipes, chapter 10. |
783 |
|
|
784 |
|
(defvar math-min-vars [(var DUMMY var-DUMMY)]) |
785 |
|
|
786 |
(defun math-min-eval (expr a) |
(defun math-min-eval (expr a) |
787 |
(if (Math-vectorp a) |
(if (Math-vectorp a) |
788 |
(let ((m -1)) |
(let ((m -1)) |
925 |
(tol (list 'float 1 (- -1 prec))) |
(tol (list 'float 1 (- -1 prec))) |
926 |
(zeps (list 'float 1 (- -5 prec))) |
(zeps (list 'float 1 (- -5 prec))) |
927 |
(e '(float 0 0)) |
(e '(float 0 0)) |
928 |
u vu xm tol1 tol2 etemp p q r xv xw) |
d u vu xm tol1 tol2 etemp p q r xv xw) |
929 |
(while (progn |
(while (progn |
930 |
(setq xm (math-mul-float '(float 5 -1) |
(setq xm (math-mul-float '(float 5 -1) |
931 |
(math-add-float a b)) |
(math-add-float a b)) |
1087 |
(list (math-add line-p xi) xi (nth 2 res)))) |
(list (math-add line-p xi) xi (nth 2 res)))) |
1088 |
|
|
1089 |
|
|
|
(defvar math-min-vars [(var DUMMY var-DUMMY)]) |
|
|
|
|
1090 |
(defun math-find-minimum (expr var guess min-widen) |
(defun math-find-minimum (expr var guess min-widen) |
1091 |
(let* ((calc-symbolic-mode nil) |
(let* ((calc-symbolic-mode nil) |
1092 |
(n 0) |
(n 0) |
1101 |
(math-dimension-error)) |
(math-dimension-error)) |
1102 |
(while (setq var (cdr var) guess (cdr guess)) |
(while (setq var (cdr var) guess (cdr guess)) |
1103 |
(or (eq (car-safe (car var)) 'var) |
(or (eq (car-safe (car var)) 'var) |
1104 |
(math-reject-arg (car vg) "*Expected a variable")) |
(math-reject-arg (car var) "*Expected a variable")) |
1105 |
(or (math-expr-contains expr (car var)) |
(or (math-expr-contains expr (car var)) |
1106 |
(math-reject-arg (car var) |
(math-reject-arg (car var) |
1107 |
"*Formula does not contain specified variable")) |
"*Formula does not contain specified variable")) |
1343 |
|
|
1344 |
|
|
1345 |
;;; Open Romberg method; "qromo" in section 4.4. |
;;; Open Romberg method; "qromo" in section 4.4. |
1346 |
|
|
1347 |
|
;; The variable math-ninteg-temp is local to math-ninteg-romberg, |
1348 |
|
;; but is used by math-ninteg-midpoint, which is used by |
1349 |
|
;; math-ninteg-romberg. |
1350 |
|
(defvar math-ninteg-temp) |
1351 |
|
|
1352 |
(defun math-ninteg-romberg (func expr lo hi mode) |
(defun math-ninteg-romberg (func expr lo hi mode) |
1353 |
(let ((curh '(float 1 0)) |
(let ((curh '(float 1 0)) |
1354 |
(h nil) |
(h nil) |
1356 |
(j 0) |
(j 0) |
1357 |
(ss nil) |
(ss nil) |
1358 |
(prec calc-internal-prec) |
(prec calc-internal-prec) |
1359 |
(integ-temp nil)) |
(math-ninteg-temp nil)) |
1360 |
(math-with-extra-prec 2 |
(math-with-extra-prec 2 |
1361 |
;; Limit on "j" loop must be 14 or less to keep "it" from overflowing. |
;; Limit on "j" loop must be 14 or less to keep "it" from overflowing. |
1362 |
(or (while (and (null ss) (<= (setq j (1+ j)) 8)) |
(or (while (and (null ss) (<= (setq j (1+ j)) 8)) |
1367 |
(if (math-lessp (math-abs (nth 1 res)) |
(if (math-lessp (math-abs (nth 1 res)) |
1368 |
(calcFunc-scf (math-abs (car res)) |
(calcFunc-scf (math-abs (car res)) |
1369 |
(- prec))) |
(- prec))) |
1370 |
(setq math-ninteg-convergence j |
(setq ss (car res))))) |
|
ss (car res))))) |
|
1371 |
(if (>= j 5) |
(if (>= j 5) |
1372 |
(setq s (cdr s) |
(setq s (cdr s) |
1373 |
h (cdr h))) |
h (cdr h))) |
1388 |
res)) |
res)) |
1389 |
|
|
1390 |
|
|
1391 |
(defun math-ninteg-midpoint (expr lo hi mode) ; uses "integ-temp" |
(defun math-ninteg-midpoint (expr lo hi mode) ; uses "math-ninteg-temp" |
1392 |
(if (eq mode 'inf) |
(if (eq mode 'inf) |
1393 |
(let ((math-infinite-mode t) temp) |
(let ((math-infinite-mode t) temp) |
1394 |
(setq temp (math-div 1 lo) |
(setq temp (math-div 1 lo) |
1395 |
lo (math-div 1 hi) |
lo (math-div 1 hi) |
1396 |
hi temp))) |
hi temp))) |
1397 |
(if integ-temp |
(if math-ninteg-temp |
1398 |
(let* ((it3 (* 3 (car integ-temp))) |
(let* ((it3 (* 3 (car math-ninteg-temp))) |
1399 |
(math-working-step-2 (* 2 (car integ-temp))) |
(math-working-step-2 (* 2 (car math-ninteg-temp))) |
1400 |
(math-working-step 0) |
(math-working-step 0) |
1401 |
(range (math-sub hi lo)) |
(range (math-sub hi lo)) |
1402 |
(del (math-div range (math-float it3))) |
(del (math-div range (math-float it3))) |
1405 |
(x (math-add lo (math-mul '(float 5 -1) del))) |
(x (math-add lo (math-mul '(float 5 -1) del))) |
1406 |
(sum '(float 0 0)) |
(sum '(float 0 0)) |
1407 |
(j 0) temp) |
(j 0) temp) |
1408 |
(while (<= (setq j (1+ j)) (car integ-temp)) |
(while (<= (setq j (1+ j)) (car math-ninteg-temp)) |
1409 |
(setq math-working-step (1+ math-working-step) |
(setq math-working-step (1+ math-working-step) |
1410 |
temp (math-ninteg-evaluate expr x mode) |
temp (math-ninteg-evaluate expr x mode) |
1411 |
math-working-step (1+ math-working-step) |
math-working-step (1+ math-working-step) |
1413 |
expr (math-add x del2) |
expr (math-add x del2) |
1414 |
mode))) |
mode))) |
1415 |
x (math-add x del3))) |
x (math-add x del3))) |
1416 |
(setq integ-temp (list it3 |
(setq math-ninteg-temp (list it3 |
1417 |
(math-add (math-div (nth 1 integ-temp) |
(math-add (math-div (nth 1 math-ninteg-temp) |
1418 |
'(float 3 0)) |
'(float 3 0)) |
1419 |
(math-mul sum del))))) |
(math-mul sum del))))) |
1420 |
(setq integ-temp (list 1 (math-mul |
(setq math-ninteg-temp (list 1 (math-mul |
1421 |
(math-sub hi lo) |
(math-sub hi lo) |
1422 |
(math-ninteg-evaluate |
(math-ninteg-evaluate |
1423 |
expr |
expr |
1424 |
(math-mul (math-add lo hi) '(float 5 -1)) |
(math-mul (math-add lo hi) '(float 5 -1)) |
1425 |
mode))))) |
mode))))) |
1426 |
(nth 1 integ-temp)) |
(nth 1 math-ninteg-temp)) |
1427 |
|
|
1428 |
|
|
1429 |
|
|
1461 |
(math-with-extra-prec 2 |
(math-with-extra-prec 2 |
1462 |
(math-general-fit expr vars coefs data 'full)))) |
(math-general-fit expr vars coefs data 'full)))) |
1463 |
|
|
1464 |
|
;; The variables math-fit-first-var, math-fit-first-coef and |
1465 |
|
;; math-fit-new-coefs are local to math-general-fit, but are used by |
1466 |
|
;; calcFunc-fitvar, calcFunc-fitparam and calcFunc-fitdummy |
1467 |
|
;; (respectively), which are used by math-general-fit. |
1468 |
|
(defvar math-fit-first-var) |
1469 |
|
(defvar math-fit-first-coef) |
1470 |
|
(defvar math-fit-new-coefs) |
1471 |
|
|
1472 |
(defun math-general-fit (expr vars coefs data mode) |
(defun math-general-fit (expr vars coefs data mode) |
1473 |
(let ((calc-simplify-mode nil) |
(let ((calc-simplify-mode nil) |
1474 |
(math-dummy-counter math-dummy-counter) |
(math-dummy-counter math-dummy-counter) |
1475 |
(math-in-fit 1) |
(math-in-fit 1) |
1476 |
(extended (eq mode 'full)) |
(extended (eq mode 'full)) |
1477 |
(first-coef math-dummy-counter) |
(math-fit-first-coef math-dummy-counter) |
1478 |
first-var |
math-fit-first-var |
1479 |
(plain-expr expr) |
(plain-expr expr) |
1480 |
orig-expr |
orig-expr |
1481 |
have-sdevs need-chisq chisq |
have-sdevs need-chisq chisq |
1483 |
(y-filter nil) |
(y-filter nil) |
1484 |
y-dummy |
y-dummy |
1485 |
(coef-filters nil) |
(coef-filters nil) |
1486 |
new-coefs |
math-fit-new-coefs |
1487 |
(xy-values nil) |
(xy-values nil) |
1488 |
(weights nil) |
(weights nil) |
1489 |
(var-YVAL nil) (var-YVALX nil) |
(var-YVAL nil) (var-YVALX nil) |
1538 |
(setq dummy (math-dummy-variable) |
(setq dummy (math-dummy-variable) |
1539 |
expr (math-expr-subst expr (car p) |
expr (math-expr-subst expr (car p) |
1540 |
(list 'calcFunc-fitparam |
(list 'calcFunc-fitparam |
1541 |
(- math-dummy-counter first-coef))))) |
(- math-dummy-counter math-fit-first-coef))))) |
1542 |
(setq first-var math-dummy-counter |
(setq math-fit-first-var math-dummy-counter |
1543 |
p vars) |
p vars) |
1544 |
(while (setq p (cdr p)) |
(while (setq p (cdr p)) |
1545 |
(or (eq (car-safe (car p)) 'var) |
(or (eq (car-safe (car p)) 'var) |
1547 |
(setq dummy (math-dummy-variable) |
(setq dummy (math-dummy-variable) |
1548 |
expr (math-expr-subst expr (car p) |
expr (math-expr-subst expr (car p) |
1549 |
(list 'calcFunc-fitvar |
(list 'calcFunc-fitvar |
1550 |
(- math-dummy-counter first-var))))) |
(- math-dummy-counter math-fit-first-var))))) |
1551 |
(if (< math-dummy-counter (+ first-var v)) |
(if (< math-dummy-counter (+ math-fit-first-var v)) |
1552 |
(setq dummy (math-dummy-variable))) ; dependent variable may be unnamed |
(setq dummy (math-dummy-variable))) ; dependent variable may be unnamed |
1553 |
(setq y-dummy dummy |
(setq y-dummy dummy |
1554 |
orig-expr expr) |
orig-expr expr) |
1607 |
(setq sigmasqr (math-add (math-sqr (nth 2 xval)) |
(setq sigmasqr (math-add (math-sqr (nth 2 xval)) |
1608 |
(or sigmasqr 0)) |
(or sigmasqr 0)) |
1609 |
xval (nth 1 xval)))) |
xval (nth 1 xval)))) |
1610 |
(set (nth 2 (aref math-dummy-vars (+ first-var j))) xval) |
(set (nth 2 (aref math-dummy-vars (+ math-fit-first-var j))) xval) |
1611 |
(setq j (1+ j))) |
(setq j (1+ j))) |
1612 |
|
|
1613 |
;; Compute Y value for this data point. |
;; Compute Y value for this data point. |
1698 |
xy-values (cdr xy-values))))) |
xy-values (cdr xy-values))))) |
1699 |
|
|
1700 |
;; Convert coefficients back into original terms. |
;; Convert coefficients back into original terms. |
1701 |
(setq new-coefs (copy-sequence beta)) |
(setq math-fit-new-coefs (copy-sequence beta)) |
1702 |
(let* ((bp new-coefs) |
(let* ((bp math-fit-new-coefs) |
1703 |
(cp covar) |
(cp covar) |
1704 |
(sigdat 1) |
(sigdat 1) |
1705 |
(math-in-fit 3) |
(math-in-fit 3) |
1715 |
(math-sqrt (math-mul (nth (setq j (1+ j)) |
(math-sqrt (math-mul (nth (setq j (1+ j)) |
1716 |
(car (setq cp (cdr cp)))) |
(car (setq cp (cdr cp)))) |
1717 |
sigdat)))))) |
sigdat)))))) |
1718 |
(setq new-coefs (math-evaluate-expr coef-filters)) |
(setq math-fit-new-coefs (math-evaluate-expr coef-filters)) |
1719 |
(if calc-fit-to-trail |
(if calc-fit-to-trail |
1720 |
(let ((bp new-coefs) |
(let ((bp math-fit-new-coefs) |
1721 |
(cp coefs) |
(cp coefs) |
1722 |
(vec nil)) |
(vec nil)) |
1723 |
(while (setq bp (cdr bp) cp (cdr cp)) |
(while (setq bp (cdr bp) cp (cdr cp)) |
1737 |
(setq vec (cons (list 'calcFunc-fitparam n) vec) |
(setq vec (cons (list 'calcFunc-fitparam n) vec) |
1738 |
n (1- n))) |
n (1- n))) |
1739 |
vec) |
vec) |
1740 |
(append (cdr new-coefs) (cdr vars)))) |
(append (cdr math-fit-new-coefs) (cdr vars)))) |
1741 |
|
|
1742 |
;; Package the result. |
;; Package the result. |
1743 |
(math-normalize |
(math-normalize |
1761 |
(defun calcFunc-fitvar (x) |
(defun calcFunc-fitvar (x) |
1762 |
(if (>= math-in-fit 2) |
(if (>= math-in-fit 2) |
1763 |
(progn |
(progn |
1764 |
(setq x (aref math-dummy-vars (+ first-var x -1))) |
(setq x (aref math-dummy-vars (+ math-fit-first-var x -1))) |
1765 |
(or (calc-var-value (nth 2 x)) x)) |
(or (calc-var-value (nth 2 x)) x)) |
1766 |
(math-reject-arg x))) |
(math-reject-arg x))) |
1767 |
|
|
1768 |
(defun calcFunc-fitparam (x) |
(defun calcFunc-fitparam (x) |
1769 |
(if (>= math-in-fit 2) |
(if (>= math-in-fit 2) |
1770 |
(progn |
(progn |
1771 |
(setq x (aref math-dummy-vars (+ first-coef x -1))) |
(setq x (aref math-dummy-vars (+ math-fit-first-coef x -1))) |
1772 |
(or (calc-var-value (nth 2 x)) x)) |
(or (calc-var-value (nth 2 x)) x)) |
1773 |
(math-reject-arg x))) |
(math-reject-arg x))) |
1774 |
|
|
1775 |
(defun calcFunc-fitdummy (x) |
(defun calcFunc-fitdummy (x) |
1776 |
(if (= math-in-fit 3) |
(if (= math-in-fit 3) |
1777 |
(nth x new-coefs) |
(nth x math-fit-new-coefs) |
1778 |
(math-reject-arg x))) |
(math-reject-arg x))) |
1779 |
|
|
1780 |
(defun calcFunc-hasfitvars (expr) |
(defun calcFunc-hasfitvars (expr) |
1801 |
(sort (mapcar 'car vars) |
(sort (mapcar 'car vars) |
1802 |
(function (lambda (x y) (string< (nth 1 x) (nth 1 y))))))) |
(function (lambda (x y) (string< (nth 1 x) (nth 1 y))))))) |
1803 |
|
|
1804 |
|
;; The variables math-all-vars-vars (the vars for math-all-vars) and |
1805 |
|
;; math-all-vars-found are local to math-all-vars-in, but are used by |
1806 |
|
;; math-all-vars-rec which is called by math-all-vars-in. |
1807 |
|
(defvar math-all-vars-vars) |
1808 |
|
(defvar math-all-vars-found) |
1809 |
|
|
1810 |
(defun math-all-vars-in (expr) |
(defun math-all-vars-in (expr) |
1811 |
(let ((vars nil) |
(let ((math-all-vars-vars nil) |
1812 |
found) |
math-all-vars-found) |
1813 |
(math-all-vars-rec expr) |
(math-all-vars-rec expr) |
1814 |
vars)) |
math-all-vars-vars)) |
1815 |
|
|
1816 |
(defun math-all-vars-rec (expr) |
(defun math-all-vars-rec (expr) |
1817 |
(if (Math-primp expr) |
(if (Math-primp expr) |
1818 |
(if (eq (car-safe expr) 'var) |
(if (eq (car-safe expr) 'var) |
1819 |
(or (math-const-var expr) |
(or (math-const-var expr) |
1820 |
(if (setq found (assoc expr vars)) |
(if (setq math-all-vars-found (assoc expr math-all-vars-vars)) |
1821 |
(setcdr found (1+ (cdr found))) |
(setcdr math-all-vars-found (1+ (cdr math-all-vars-found))) |
1822 |
(setq vars (cons (cons expr 1) vars))))) |
(setq math-all-vars-vars (cons (cons expr 1) math-all-vars-vars))))) |
1823 |
(while (setq expr (cdr expr)) |
(while (setq expr (cdr expr)) |
1824 |
(math-all-vars-rec (car expr))))) |
(math-all-vars-rec (car expr))))) |
1825 |
|
|
1826 |
|
(provide 'calcalg3) |
1827 |
|
|
1828 |
;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6 |
;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6 |
1829 |
;;; calcalg3.el ends here |
;;; calcalg3.el ends here |