Commit 03cc1aba authored by Jay Belanger's avatar Jay Belanger
Browse files

(calc-curve-nvars, calc-curve-varnames, calc-curve-model)

(calc-curve-coefnames): New variable.
(calc-curve-fit, calc-get-fit-variables): Replace variables nvars,
varnames, model and coefnames by declared variables.

(math-root-widen): New variable.
(math-search-root, math-find-root): Replace variable root-widen by
declared variable.

(var-DUMMY): Declare it.
(math-root-vars, math-min-vars): Move the declarations to earlier in
the file.

(math-brent-min): Make d a local variable.

(math-find-minimum): Replace non-existent variable.

(math-ninteg-romberg): Remove unnecessary variable.

(math-ninteg-temp): New variable.
(math-ninteg-romberg, math-ninteg-midpoint): Replace variable
integ-temp by declared variable.

(math-fit-first-var, math-fit-first-coef, math-fit-new-coefs): New variables.
(math-general-fit): Replace variables first-var, first-coef and
new-coefs by declared variables.
(calcFunc-fitvar): Replace variable first-var by declared variable.
(calcFunc-fitparam): Replace variable first-coef by declared variable.
(calcFunc-fitdummy): Replace variable new-coefs by declared variable.

(math-all-vars-vars, math-all-vars-found): New variables.
(math-all-vars-in, math-all-vars-rec): Replace variables vars and
found by declared variable math-all-vars-vars.
parent a6cecab9
......@@ -3,8 +3,7 @@
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
;; Author: David Gillespie <>
;; Maintainers: D. Goel <>
;; Colin Walters <>
;; Maintainer: Jay Belanger <>
;; This file is part of GNU Emacs.
......@@ -99,8 +98,15 @@
(calc-enter-result 1 "poli" (list 'calcFunc-polint data
(calc-top 1)))))))
;; The variables calc-curve-nvars, calc-curve-varnames, calc-curve-model and calc-curve-coefnames are local to calc-curve-fit, but are
;; used by calc-get-fit-variables which is called by calc-curve-fit.
(defvar calc-curve-nvars)
(defvar calc-curve-varnames)
(defvar calc-curve-model)
(defvar calc-curve-coefnames)
(defun calc-curve-fit (arg &optional model coefnames varnames)
(defun calc-curve-fit (arg &optional calc-curve-model
calc-curve-coefnames calc-curve-varnames)
(interactive "P")
(setq calc-aborted-prefix nil)
......@@ -108,7 +114,7 @@
(if (calc-is-hyperbolic) 'calcFunc-efit
key (which 0)
n nvars temp data
n calc-curve-nvars temp data
(homog nil)
(msgs '( "(Press ? for help)"
"1 = linear or multilinear"
......@@ -120,7 +126,7 @@
"g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
"h prefix = homogeneous model (no constant term)"
"' = alg entry, $ = stack, u = Model1, U = Model2")))
(while (not model)
(while (not calc-curve-model)
(message "Fit to model: %s:%s"
(nth which msgs)
(if homog " h" ""))
......@@ -150,44 +156,50 @@
(t (error "Bad prefix argument")))
(or (math-matrixp data) (not (cdr (cdr data)))
(error "Data matrix is not a matrix!"))
(setq nvars (- (length data) 2)
coefnames nil
varnames nil)
(setq calc-curve-nvars (- (length data) 2)
calc-curve-coefnames nil
calc-curve-varnames nil)
((= key ?1) ; linear or multilinear
(calc-get-fit-variables nvars (1+ nvars) (and homog 0))
(setq model (math-mul coefnames
(cons 'vec (cons 1 (cdr varnames))))))
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 0))
(setq calc-curve-model (math-mul calc-curve-coefnames
(cons 'vec (cons 1 (cdr calc-curve-varnames))))))
((and (>= key ?2) (<= key ?9)) ; polynomial
(calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
(setq model (math-build-polynomial-expr (cdr coefnames)
(nth 1 varnames))))
(setq calc-curve-model
(math-build-polynomial-expr (cdr calc-curve-coefnames)
(nth 1 calc-curve-varnames))))
((= key ?i) ; exact polynomial
(calc-get-fit-variables 1 (1- (length (nth 1 data)))
(and homog 0))
(setq model (math-build-polynomial-expr (cdr coefnames)
(nth 1 varnames))))
(setq calc-curve-model
(math-build-polynomial-expr (cdr calc-curve-coefnames)
(nth 1 calc-curve-varnames))))
((= key ?p) ; power law
(calc-get-fit-variables nvars (1+ nvars) (and homog 1))
(setq model (math-mul (nth 1 coefnames)
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 1))
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
'(var mul var-mul)
'(var pow var-pow)
(cons 'vec (cdr (cdr coefnames))))))))
(cons 'vec (cdr (cdr calc-curve-coefnames))))))))
((= key ?^) ; exponential law
(calc-get-fit-variables nvars (1+ nvars) (and homog 1))
(setq model (math-mul (nth 1 coefnames)
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 1))
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
'(var mul var-mul)
'(var pow var-pow)
(cons 'vec (cdr (cdr coefnames)))
(cons 'vec (cdr (cdr calc-curve-coefnames)))
((memq key '(?e ?E))
(calc-get-fit-variables nvars (1+ nvars) (and homog 1))
(setq model (math-mul (nth 1 coefnames)
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 1))
(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
'(var mul var-mul)
......@@ -198,45 +210,50 @@
(^ 10 (var a var-a))))
'(var mul var-mul)
(cons 'vec (cdr (cdr coefnames)))
(cons 'vec (cdr (cdr calc-curve-coefnames)))
((memq key '(?x ?X))
(calc-get-fit-variables nvars (1+ nvars) (and homog 0))
(setq model (math-mul coefnames
(cons 'vec (cons 1 (cdr varnames)))))
(setq model (if (eq key ?x)
(list 'calcFunc-exp model)
(list '^ 10 model))))
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 0))
(setq calc-curve-model (math-mul calc-curve-coefnames
(cons 'vec (cons 1 (cdr calc-curve-varnames)))))
(setq calc-curve-model (if (eq key ?x)
(list 'calcFunc-exp calc-curve-model)
(list '^ 10 calc-curve-model))))
((memq key '(?l ?L))
(calc-get-fit-variables nvars (1+ nvars) (and homog 0))
(setq model (math-mul coefnames
(calc-get-fit-variables calc-curve-nvars
(1+ calc-curve-nvars) (and homog 0))
(setq calc-curve-model (math-mul calc-curve-coefnames
(cons 'vec
(cons 1 (cdr (calcFunc-map
(if (eq key ?l)
'(var ln var-ln)
'(var log10
((= key ?q)
(calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0))
(let ((c coefnames)
(v varnames))
(setq model (nth 1 c))
(calc-get-fit-variables calc-curve-nvars
(1+ (* 2 calc-curve-nvars)) (and homog 0))
(let ((c calc-curve-coefnames)
(v calc-curve-varnames))
(setq calc-curve-model (nth 1 c))
(while (setq v (cdr v) c (cdr (cdr c)))
(setq model (math-add
(setq calc-curve-model (math-add
(list '*
(car c)
(list '^
(list '- (car v) (nth 1 c))
((= key ?g)
(setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
varnames '(vec (var XFit var-XFit))
coefnames '(vec (var AFit var-AFit)
(setq calc-curve-model
(math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
calc-curve-varnames '(vec (var XFit var-XFit))
calc-curve-coefnames '(vec (var AFit var-AFit)
(var BFit var-BFit)
(var CFit var-CFit)))
(calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))
(calc-get-fit-variables 1 (1- (length calc-curve-coefnames))
(and homog 1)))
((memq key '(?\$ ?\' ?u ?U))
(let* ((defvars nil)
(record-entry nil))
......@@ -244,74 +261,78 @@
(let* ((calc-dollar-values calc-arg-values)
(calc-dollar-used 0)
(calc-hashes-used 0))
(setq model (calc-do-alg-entry "" "Model formula: "))
(if (/= (length model) 1)
(setq calc-curve-model (calc-do-alg-entry "" "Model formula: "))
(if (/= (length calc-curve-model) 1)
(error "Bad format"))
(setq model (car model)
(setq calc-curve-model (car calc-curve-model)
record-entry t)
(if (> calc-dollar-used 0)
(setq coefnames
(setq calc-curve-coefnames
(cons 'vec
(nthcdr (- (length calc-arg-values)
(reverse calc-arg-values))))
(if (> calc-hashes-used 0)
(setq coefnames
(setq calc-curve-coefnames
(cons 'vec (calc-invent-args
(setq model (cond ((eq key ?u)
(setq calc-curve-model (cond ((eq key ?u)
(calc-var-value 'var-Model1))
((eq key ?U)
(calc-var-value 'var-Model2))
(t (calc-top 1))))
(or model (error "User model not yet defined"))
(if (math-vectorp model)
(if (and (memq (length model) '(3 4))
(not (math-objvecp (nth 1 model)))
(math-vectorp (nth 2 model))
(or (null (nth 3 model))
(math-vectorp (nth 3 model))))
(setq varnames (nth 2 model)
coefnames (or (nth 3 model)
(or calc-curve-model (error "User model not yet defined"))
(if (math-vectorp calc-curve-model)
(if (and (memq (length calc-curve-model) '(3 4))
(not (math-objvecp (nth 1 calc-curve-model)))
(math-vectorp (nth 2 calc-curve-model))
(or (null (nth 3 calc-curve-model))
(math-vectorp (nth 3 calc-curve-model))))
(setq calc-curve-varnames (nth 2 calc-curve-model)
(or (nth 3 calc-curve-model)
(cons 'vec
model varnames)))
model (nth 1 model))
calc-curve-model calc-curve-varnames)))
calc-curve-model (nth 1 calc-curve-model))
(error "Incorrect model specifier")))))
(or varnames
(let ((with-y (eq (car-safe model) 'calcFunc-eq)))
(if coefnames
(calc-get-fit-variables (if with-y (1+ nvars) nvars)
(1- (length coefnames))
(or calc-curve-varnames
(let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq)))
(if calc-curve-coefnames
(if with-y (1+ calc-curve-nvars) calc-curve-nvars)
(1- (length calc-curve-coefnames))
model coefnames)
calc-curve-model calc-curve-coefnames)
nil with-y)
(let* ((coefs (math-all-vars-but model nil))
(let* ((coefs (math-all-vars-but calc-curve-model nil))
(vars nil)
(n (- (length coefs) nvars (if with-y 2 1)))
(n (- (length coefs) calc-curve-nvars (if with-y 2 1)))
(if (< n 0)
(error "Not enough variables in model"))
(setq p (nthcdr n coefs))
(setq vars (cdr p))
(setcdr p nil)
(calc-get-fit-variables (if with-y (1+ nvars) nvars)
(if with-y (1+ calc-curve-nvars) calc-curve-nvars)
(length coefs)
vars coefs with-y)))))
(if record-entry
(calc-record (list 'vec model varnames coefnames)
(calc-record (list 'vec calc-curve-model
calc-curve-varnames calc-curve-coefnames)
(t (beep))))
(let ((calc-fit-to-trail t))
(calc-enter-result n (substring (symbol-name func) 9)
(list func model
(if (= (length varnames) 2)
(nth 1 varnames)
(if (= (length coefnames) 2)
(nth 1 coefnames)
(list func calc-curve-model
(if (= (length calc-curve-varnames) 2)
(nth 1 calc-curve-varnames)
(if (= (length calc-curve-coefnames) 2)
(nth 1 calc-curve-coefnames)
(if (consp calc-fit-to-trail)
(calc-record (calc-normalize calc-fit-to-trail) "parm"))))))
......@@ -340,7 +361,7 @@
(calc-invent-variables num but t base))))
(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
(or (= nv (if with-y (1+ nvars) nvars))
(or (= nv (if with-y (1+ calc-curve-nvars) calc-curve-nvars))
(error "Wrong number of data vectors for this type of model"))
(if (integerp defv)
(setq homog defv
......@@ -388,12 +409,12 @@
(error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
(if homog
(setq coefs (cons 'vec (cons homog (cdr coefs)))))
(if varnames
(setq model (math-multi-subst model (cdr varnames) (cdr vars))))
(if coefnames
(setq model (math-multi-subst model (cdr coefnames) (cdr coefs))))
(setq varnames vars
coefnames coefs)))
(if calc-curve-varnames
(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-varnames) (cdr vars))))
(if calc-curve-coefnames
(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-coefnames) (cdr coefs))))
(setq calc-curve-varnames vars
calc-curve-coefnames coefs)))
......@@ -401,6 +422,9 @@
;;; The following algorithms are from Numerical Recipes chapter 9.
;;; "rtnewt" with safety kludges
(defvar var-DUMMY)
(defun math-newton-root (expr deriv guess orig-guess limit)
(math-working "newton" guess)
(let* ((var-DUMMY guess)
......@@ -494,14 +518,20 @@
low vlow high vhigh))))))
;;; Search for a root in an interval with no overt zero crossing.
;; The variable math-root-widen is local to math-find-root, but
;; is used by math-search-root, which is called (directly and
;; indirectly) by math-find-root.
(defvar math-root-widen)
(defun math-search-root (expr deriv low vlow high vhigh)
(let (found)
(if root-widen
(if math-root-widen
(let ((iters 0)
(iterlim (if (eq root-widen 'point)
(iterlim (if (eq math-root-widen 'point)
(+ calc-internal-prec 10)
(factor (if (eq root-widen 'point)
(factor (if (eq math-root-widen 'point)
'(float 9 0)
'(float 16 -1)))
(prev nil) vprev waslow
......@@ -600,6 +630,9 @@
(list 'vec mid vmid)))
;;; "mnewt"
(defvar math-root-vars [(var DUMMY var-DUMMY)])
(defun math-newton-multi (expr jacob n guess orig-guess limit)
(let ((m -1)
(p guess)
......@@ -624,9 +657,8 @@
(math-reject-arg nil "*Newton's method failed to converge"))
(list 'vec next expr-val))))
(defvar math-root-vars [(var DUMMY var-DUMMY)])
(defun math-find-root (expr var guess root-widen)
(defun math-find-root (expr var guess math-root-widen)
(if (eq (car-safe expr) 'vec)
(let ((n (1- (length expr)))
(calc-symbolic-mode nil)
......@@ -710,7 +742,7 @@
var-DUMMY guess
vlow (math-evaluate-expr expr)
vhigh vlow
root-widen 'point)
math-root-widen 'point)
(if (eq (car guess) 'intv)
(or (math-constp guess) (math-reject-arg guess 'constp))
......@@ -752,6 +784,8 @@
;;; The following algorithms come from Numerical Recipes, chapter 10.
(defvar math-min-vars [(var DUMMY var-DUMMY)])
(defun math-min-eval (expr a)
(if (Math-vectorp a)
(let ((m -1))
......@@ -894,7 +928,7 @@
(tol (list 'float 1 (- -1 prec)))
(zeps (list 'float 1 (- -5 prec)))
(e '(float 0 0))
u vu xm tol1 tol2 etemp p q r xv xw)
d u vu xm tol1 tol2 etemp p q r xv xw)
(while (progn
(setq xm (math-mul-float '(float 5 -1)
(math-add-float a b))
......@@ -1056,8 +1090,6 @@
(list (math-add line-p xi) xi (nth 2 res))))
(defvar math-min-vars [(var DUMMY var-DUMMY)])
(defun math-find-minimum (expr var guess min-widen)
(let* ((calc-symbolic-mode nil)
(n 0)
......@@ -1072,7 +1104,7 @@
(while (setq var (cdr var) guess (cdr guess))
(or (eq (car-safe (car var)) 'var)
(math-reject-arg (car vg) "*Expected a variable"))
(math-reject-arg (car var) "*Expected a variable"))
(or (math-expr-contains expr (car var))
(math-reject-arg (car var)
"*Formula does not contain specified variable"))
......@@ -1314,6 +1346,12 @@
;;; Open Romberg method; "qromo" in section 4.4.
;; The variable math-ninteg-temp is local to math-ninteg-romberg,
;; but is used by math-ninteg-midpoint, which is used by
;; math-ninteg-romberg.
(defvar math-ninteg-temp)
(defun math-ninteg-romberg (func expr lo hi mode)
(let ((curh '(float 1 0))
(h nil)
......@@ -1321,7 +1359,7 @@
(j 0)
(ss nil)
(prec calc-internal-prec)
(integ-temp nil))
(math-ninteg-temp nil))
(math-with-extra-prec 2
;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.
(or (while (and (null ss) (<= (setq j (1+ j)) 8))
......@@ -1332,8 +1370,7 @@
(if (math-lessp (math-abs (nth 1 res))
(calcFunc-scf (math-abs (car res))
(- prec)))
(setq math-ninteg-convergence j
ss (car res)))))
(setq ss (car res)))))
(if (>= j 5)
(setq s (cdr s)
h (cdr h)))
......@@ -1354,15 +1391,15 @@
(defun math-ninteg-midpoint (expr lo hi mode) ; uses "integ-temp"
(defun math-ninteg-midpoint (expr lo hi mode) ; uses "math-ninteg-temp"
(if (eq mode 'inf)
(let ((math-infinite-mode t) temp)
(setq temp (math-div 1 lo)
lo (math-div 1 hi)
hi temp)))
(if integ-temp
(let* ((it3 (* 3 (car integ-temp)))
(math-working-step-2 (* 2 (car integ-temp)))
(if math-ninteg-temp
(let* ((it3 (* 3 (car math-ninteg-temp)))
(math-working-step-2 (* 2 (car math-ninteg-temp)))
(math-working-step 0)
(range (math-sub hi lo))
(del (math-div range (math-float it3)))
......@@ -1371,7 +1408,7 @@
(x (math-add lo (math-mul '(float 5 -1) del)))
(sum '(float 0 0))
(j 0) temp)
(while (<= (setq j (1+ j)) (car integ-temp))
(while (<= (setq j (1+ j)) (car math-ninteg-temp))
(setq math-working-step (1+ math-working-step)
temp (math-ninteg-evaluate expr x mode)
math-working-step (1+ math-working-step)
......@@ -1379,17 +1416,17 @@
expr (math-add x del2)
x (math-add x del3)))
(setq integ-temp (list it3
(math-add (math-div (nth 1 integ-temp)
(setq math-ninteg-temp (list it3
(math-add (math-div (nth 1 math-ninteg-temp)
'(float 3 0))
(math-mul sum del)))))
(setq integ-temp (list 1 (math-mul
(setq math-ninteg-temp (list 1 (math-mul
(math-sub hi lo)
(math-mul (math-add lo hi) '(float 5 -1))
(nth 1 integ-temp))
(nth 1 math-ninteg-temp))
......@@ -1427,13 +1464,21 @@
(math-with-extra-prec 2
(math-general-fit expr vars coefs data 'full))))
;; The variables math-fit-first-var, math-fit-first-coef and
;; math-fit-new-coefs are local to math-general-fit, but are used by
;; calcFunc-fitvar, calcFunc-fitparam and calcFunc-fitdummy
;; (respectively), which are used by math-general-fit.
(defvar math-fit-first-var)
(defvar math-fit-first-coef)
(defvar math-fit-new-coefs)
(defun math-general-fit (expr vars coefs data mode)
(let ((calc-simplify-mode nil)
(math-dummy-counter math-dummy-counter)
(math-in-fit 1)
(extended (eq mode 'full))
(first-coef math-dummy-counter)
(math-fit-first-coef math-dummy-counter)
(plain-expr expr)
have-sdevs need-chisq chisq
......@@ -1441,7 +1486,7 @@
(y-filter nil)
(coef-filters nil)
(xy-values nil)
(weights nil)
(var-YVAL nil) (var-YVALX nil)
......@@ -1496,8 +1541,8 @@
(setq dummy (math-dummy-variable)
expr (math-expr-subst expr (car p)
(list 'calcFunc-fitparam
(- math-dummy-counter first-coef)))))
(setq first-var math-dummy-counter
(- math-dummy-counter math-fit-first-coef)))))
(setq math-fit-first-var math-dummy-counter
p vars)
(while (setq p (cdr p))
(or (eq (car-safe (car p)) 'var)
......@@ -1505,8 +1550,8 @@
(setq dummy (math-dummy-variable)
expr (math-expr-subst expr (car p)
(list 'calcFunc-fitvar
(- math-dummy-counter first-var)))))
(if (< math-dummy-counter (+ first-var v))
(- math-dummy-counter math-fit-first-var)))))
(if (< math-dummy-counter (+ math-fit-first-var v))
(setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
(setq y-dummy dummy
orig-expr expr)
......@@ -1565,7 +1610,7 @@
(setq sigmasqr (math-add (math-sqr (nth 2 xval))
(or sigmasqr 0))
xval (nth 1 xval))))
(set (nth 2 (aref math-dummy-vars (+ first-var j))) xval)
(set (nth 2 (aref math-dummy-vars (+ math-fit-first-var j))) xval)
(setq j (1+ j)))
;; Compute Y value for this data point.
......@@ -1656,8 +1701,8 @@
xy-values (cdr xy-values)))))
;; Convert coefficients back into original terms.
(setq new-coefs (copy-sequence beta))
(let* ((bp new-coefs)
(setq math-fit-new-coefs (copy-sequence beta))
(let* ((bp math-fit-new-coefs)
(cp covar)
(sigdat 1)
(math-in-fit 3)
......@@ -1673,9 +1718,9 @@
(math-sqrt (math-mul (nth (setq j (1+ j))
(car (setq cp (cdr cp))))
(setq new-coefs (math-evaluate-expr coef-filters))
(setq math-fit-new-coefs (math-evaluate-expr coef-filters))
(if calc-fit-to-trail
(let ((bp new-coefs)
(let ((bp math-fit-new-coefs)
(cp coefs)
(vec nil))
(while (setq bp (cdr bp) cp (cdr cp))
......@@ -1695,7 +1740,7 @@
(setq vec (cons (list 'calcFunc-fitparam n) vec)
n (1- n)))
(append (cdr new-coefs) (cdr vars))))
(append (cdr math-fit-new-coefs) (cdr vars))))
;; Package the result.
......@@ -1719,20 +1764,20 @@
(defun calcFunc-fitvar (x)
(if (>= math-in-fit 2)
(setq x (aref math-dummy-vars (+ first-var x -1)))
(setq x (aref math-dummy-vars (+ math-fit-first-var x -1)))
(or (calc-var-value (nth 2 x)) x))
(math-reject-arg x)))
(defun calcFunc-fitparam (x)
(if (>= math-in-fit 2)
(setq x (aref math-dummy-vars (+ first-coef x -1)))
(setq x (aref math-dummy-vars (+ math-fit-first-coef x -1)))
(or (calc-var-value (nth 2 x)) x))
(math-reject-arg x)))
(defun calcFunc-fitdummy (x)
(if (= math-in-fit 3)
(nth x new-coefs)
(nth x math-fit-new-coefs)