-;; Calculator for GNU Emacs, part II [calc-alg-3.el]
-;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
-;; Written by Dave Gillespie, daveg@synaptics.com.
+;;; calcalg3.el --- more algebraic functions for Calc
+
+;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
+;; 2005, 2006, 2007 Free Software Foundation, Inc.
+
+;; Author: David Gillespie <daveg@synaptics.com>
+;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
;; This file is part of GNU Emacs.
+;; GNU Emacs is free software; you can redistribute it and/or modify
+;; it under the terms of the GNU General Public License as published by
+;; the Free Software Foundation; either version 3, or (at your option)
+;; any later version.
+
;; GNU Emacs is distributed in the hope that it will be useful,
-;; but WITHOUT ANY WARRANTY. No author or distributor
-;; accepts responsibility to anyone for the consequences of using it
-;; or for whether it serves any particular purpose or works at all,
-;; unless he says so in writing. Refer to the GNU Emacs General Public
-;; License for full details.
+;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+;; GNU General Public License for more details.
-;; Everyone is granted permission to copy, modify and redistribute
-;; GNU Emacs, but only under the conditions described in the
-;; GNU Emacs General Public License. A copy of this license is
-;; supposed to have been given to you along with GNU Emacs so you
-;; can know your rights and responsibilities. It should be in a
-;; file named COPYING. Among other things, the copyright notice
-;; and this notice must be preserved on all copies.
+;; You should have received a copy of the GNU General Public License
+;; along with GNU Emacs; see the file COPYING. If not, write to the
+;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+;; Boston, MA 02110-1301, USA.
+;;; Commentary:
+;;; Code:
;; This file is autoloaded from calc-ext.el.
-(require 'calc-ext)
+(require 'calc-ext)
(require 'calc-macs)
-(defun calc-Need-calc-alg-3 () nil)
+;; Declare functions which are defined elsewhere.
+(declare-function calc-fit-s-shaped-logistic-curve "calc-nlfit" (arg))
+(declare-function calc-fit-bell-shaped-logistic-curve "calc-nlfit" (arg))
+(declare-function calc-fit-hubbert-linear-curve "calc-nlfit" (&optional sdv))
+(declare-function calc-graph-add-curve "calc-graph" (xdata ydata &optional zdata))
+(declare-function calc-graph-lookup "calc-graph" (thing))
+(declare-function calc-graph-set-styles "calc-graph" (lines points &optional yerr))
+(declare-function math-min-list "calc-arith" (a b))
+(declare-function math-max-list "calc-arith" (a b))
+
+(defun math-map-binop (binop args1 args2)
+ "Apply BINOP to the elements of the lists ARGS1 and ARGS2"
+ (if args1
+ (cons
+ (funcall binop (car args1) (car args2))
+ (funcall 'math-map-binop binop (cdr args1) (cdr args2)))))
(defun calc-find-root (var)
(interactive "sVariable(s) to solve for: ")
(calc-enter-result 1 "root" (list func
(calc-top-n 2)
var
- (calc-top-n 1)))))))
-)
+ (calc-top-n 1))))))))
(defun calc-find-minimum (var)
(interactive "sVariable(s) to minimize over: ")
(calc-enter-result 1 tag (list func
(calc-top-n 2)
var
- (calc-top-n 1)))))))
-)
+ (calc-top-n 1))))))))
(defun calc-find-maximum (var)
(interactive "sVariable to maximize over: ")
(calc-invert-func)
- (calc-find-minimum var)
-)
+ (calc-find-minimum var))
(defun calc-poly-interp (arg)
(if (calc-is-hyperbolic)
(calc-enter-result 1 "rati" (list 'calcFunc-ratint data (calc-top 1)))
(calc-enter-result 1 "poli" (list 'calcFunc-polint data
- (calc-top 1))))))
-)
+ (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)
+(defvar calc-curve-fit-history nil
+ "History for calc-curve-fit.")
+
+(defun calc-curve-fit (arg &optional calc-curve-model
+ calc-curve-coefnames calc-curve-varnames)
(interactive "P")
(calc-slow-wrapper
(setq calc-aborted-prefix nil)
(if (calc-is-hyperbolic) 'calcFunc-efit
'calcFunc-fit)))
key (which 0)
- n nvars temp data
+ (nonlinear nil)
+ (plot nil)
+ n calc-curve-nvars temp data
(homog nil)
(msgs '( "(Press ? for help)"
"1 = linear or multilinear"
"E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
"q = a + b (x-c)^2"
"g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
+ "s = a/(1 + exp(b (x - c)))"
+ "b = a exp(b (x - c))/(1 + exp(b (x - c)))^2"
+ "o = (y/x) = a (1 - x/b)"
"h prefix = homogeneous model (no constant term)"
+ "P prefix = plot result"
"' = alg entry, $ = stack, u = Model1, U = Model2")))
- (while (not model)
- (message "Fit to model: %s:%s"
- (nth which msgs)
- (if homog " h" ""))
+ (while (not calc-curve-model)
+ (message
+ "Fit to model: %s:%s%s"
+ (nth which msgs)
+ (if plot "P" " ")
+ (if homog "h" ""))
(setq key (read-char))
(cond ((= key ?\C-g)
(keyboard-quit))
(setq which (% (1+ which) (length msgs))))
((memq key '(?h ?H))
(setq homog (not homog)))
+ ((= key ?P)
+ (if plot
+ (setq plot nil)
+ (let ((data (calc-top 1)))
+ (if (or
+ (calc-is-hyperbolic)
+ (calc-is-inverse)
+ (not (= (length data) 3)))
+ (setq plot "Can't plot")
+ (setq plot data)))))
((progn
(if (eq key ?\$)
(setq n 1)
(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)
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)
- (calcFunc-reduce
- '(var mul var-mul)
- (calcFunc-map
- '(var pow var-pow)
- varnames
- (cons 'vec (cdr (cdr 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)
+ (calcFunc-reduce
+ '(var mul var-mul)
+ (calcFunc-map
+ '(var pow var-pow)
+ calc-curve-varnames
+ (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)
- (calcFunc-reduce
- '(var mul var-mul)
- (calcFunc-map
- '(var pow var-pow)
- (cons 'vec (cdr (cdr coefnames)))
- varnames)))))
+ (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)
+ (calcFunc-reduce
+ '(var mul var-mul)
+ (calcFunc-map
+ '(var pow var-pow)
+ (cons 'vec (cdr (cdr calc-curve-coefnames)))
+ calc-curve-varnames)))))
+ ((= key ?s)
+ (setq nonlinear t)
+ (setq calc-curve-model t)
+ (require 'calc-nlfit)
+ (calc-fit-s-shaped-logistic-curve func))
+ ((= key ?b)
+ (setq nonlinear t)
+ (setq calc-curve-model t)
+ (require 'calc-nlfit)
+ (calc-fit-bell-shaped-logistic-curve func))
+ ((= key ?o)
+ (setq nonlinear t)
+ (setq calc-curve-model t)
+ (require 'calc-nlfit)
+ (if (and plot (not (stringp plot)))
+ (setq plot
+ (list 'vec
+ (nth 1 plot)
+ (cons
+ 'vec
+ (math-map-binop 'calcFunc-div
+ (cdr (nth 2 plot))
+ (cdr (nth 1 plot)))))))
+ (calc-fit-hubbert-linear-curve func))
((memq key '(?e ?E))
- (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
- (setq model (math-mul (nth 1 coefnames)
- (calcFunc-reduce
- '(var mul var-mul)
- (calcFunc-map
- (if (eq key ?e)
- '(var exp var-exp)
- '(calcFunc-lambda
- (var a var-a)
- (^ 10 (var a var-a))))
- (calcFunc-map
- '(var mul var-mul)
- (cons 'vec (cdr (cdr coefnames)))
- varnames))))))
+ (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)
+ (calcFunc-reduce
+ '(var mul var-mul)
+ (calcFunc-map
+ (if (eq key ?e)
+ '(var exp var-exp)
+ '(calcFunc-lambda
+ (var a var-a)
+ (^ 10 (var a var-a))))
+ (calcFunc-map
+ '(var mul var-mul)
+ (cons 'vec (cdr (cdr calc-curve-coefnames)))
+ calc-curve-varnames))))))
((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
- (cons 'vec
- (cons 1 (cdr (calcFunc-map
- (if (eq key ?l)
- '(var ln var-ln)
- '(var log10
- var-log10))
- 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 (calcFunc-map
+ (if (eq key ?l)
+ '(var ln var-ln)
+ '(var log10
+ var-log10))
+ calc-curve-varnames)))))))
((= 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
- model
+ (setq calc-curve-model (math-add
+ calc-curve-model
(list '*
(car c)
(list '^
(list '- (car v) (nth 1 c))
2)))))))
((= 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)
- (var BFit var-BFit)
- (var CFit var-CFit)))
- (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))
+ (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 calc-curve-coefnames))
+ (and homog 1)))
((memq key '(?\$ ?\' ?u ?U))
(let* ((defvars nil)
(record-entry nil))
(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: "
+ nil 'calc-curve-fit-history))
+ (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)
calc-dollar-used)
(reverse calc-arg-values))))
(if (> calc-hashes-used 0)
- (setq coefnames
+ (setq calc-curve-coefnames
(cons 'vec (calc-invent-args
calc-hashes-used))))))
(progn
- (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)
- (cons 'vec
- (math-all-vars-but
- model varnames)))
- model (nth 1 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)
+ calc-curve-coefnames
+ (or (nth 3 calc-curve-model)
+ (cons 'vec
+ (math-all-vars-but
+ 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))
- (math-all-vars-but
- model coefnames)
- nil with-y)
- (let* ((coefs (math-all-vars-but model nil))
+ (or calc-curve-varnames
+ (let ((with-y
+ (eq (car-safe calc-curve-model) 'calcFunc-eq)))
+ (if calc-curve-coefnames
+ (calc-get-fit-variables
+ (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
+ (1- (length calc-curve-coefnames))
+ (math-all-vars-but
+ calc-curve-model calc-curve-coefnames)
+ nil with-y)
+ (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)))
p)
(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)
- (length coefs)
- vars coefs with-y)))))
+ (calc-get-fit-variables
+ (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)
"modl"))))
(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)
- varnames)
- (if (= (length coefnames) 2)
- (nth 1 coefnames)
- coefnames)
- data))
- (if (consp calc-fit-to-trail)
- (calc-record (calc-normalize calc-fit-to-trail) "parm")))))
-)
+ (unless nonlinear
+ (let ((calc-fit-to-trail t))
+ (calc-enter-result n (substring (symbol-name func) 9)
+ (list func calc-curve-model
+ (if (= (length calc-curve-varnames) 2)
+ (nth 1 calc-curve-varnames)
+ calc-curve-varnames)
+ (if (= (length calc-curve-coefnames) 2)
+ (nth 1 calc-curve-coefnames)
+ calc-curve-coefnames)
+ data))
+ (if (consp calc-fit-to-trail)
+ (calc-record (calc-normalize calc-fit-to-trail) "parm"))))
+ (when plot
+ (if (stringp plot)
+ (message "%s" plot)
+ (let ((calc-graph-no-auto-view t))
+ (calc-graph-delete t)
+ (calc-graph-add-curve
+ (calc-graph-lookup (nth 1 plot))
+ (calc-graph-lookup (nth 2 plot)))
+ (unless (math-contains-sdev-p (nth 2 data))
+ (calc-graph-set-styles nil nil)
+ (calc-graph-point-style nil))
+ (setq plot (cdr (nth 1 plot)))
+ (setq plot
+ (list 'intv
+ 3
+ (math-sub
+ (math-min-list (car plot) (cdr plot))
+ '(float 5 -1))
+ (math-add
+ '(float 5 -1)
+ (math-max-list (car plot) (cdr plot)))))
+ (calc-graph-add-curve (calc-graph-lookup plot)
+ (calc-graph-lookup (calc-top-n 1)))
+ (calc-graph-plot nil)))))))
(defun calc-invent-independent-variables (n &optional but)
- (calc-invent-variables n but '(x y z t) "x")
-)
+ (calc-invent-variables n but '(x y z t) "x"))
(defun calc-invent-parameter-variables (n &optional but)
- (calc-invent-variables n but '(a b c d) "a")
-)
+ (calc-invent-variables n but '(a b c d) "a"))
(defun calc-invent-variables (num but names base)
(let ((vars nil)
(while (and (> n 0) names)
(setq var (math-build-var-name (if (consp names)
(car names)
- (concat base (setq nn (1+ nn))))))
+ (concat base (int-to-string
+ (setq nn (1+ nn)))))))
(or (math-expr-contains (cons 'vec but) var)
(setq vars (cons var vars)
n (1- n)))
(or (symbolp names) (setq names (cdr names))))
(if (= n 0)
(nreverse vars)
- (calc-invent-variables num but t base)))
-)
+ (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
(setq defv (calc-invent-independent-variables nv)))
(or defc
(setq defc (calc-invent-parameter-variables nc defv)))
- (let ((vars (read-string (format "Fitting variables: (default %s; %s) "
+ (let ((vars (read-string (format "Fitting variables (default %s; %s): "
(mapconcat 'symbol-name
(mapcar (function (lambda (v)
(nth 1 v)))
(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)))
;;; 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)
limit)
(math-newton-root expr deriv next orig-guess limit)
(math-reject-arg next "*Newton's method failed to converge"))))
- (math-reject-arg next "*Newton's method encountered a singularity")))
-)
+ (math-reject-arg next "*Newton's method encountered a singularity"))))
;;; Inspired by "rtsafe"
(defun math-newton-search-root (expr deriv guess vguess ostep oostep
(and (Math-negp vlow) (Math-negp vhigh)))
(math-search-root expr deriv low vlow high vhigh)
(math-newton-search-root expr deriv nil nil nil ostep
- low vlow high vhigh)))))
-)
+ 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)
20))
- (factor (if (eq root-widen 'point)
+ (factor (if (eq math-root-widen 'point)
'(float 9 0)
'(float 16 -1)))
(prev nil) vprev waslow
low vlow high vhigh)
(math-bisect-root expr low vlow high vhigh))))
(math-reject-arg (list 'intv 3 low high)
- "*Unable to find a sign change in this interval")))
-)
+ "*Unable to find a sign change in this interval"))))
;;; "rtbis" (but we should be using Brent's method)
(defun math-bisect-root (expr low vlow high vhigh)
vhigh vmid)
(setq low mid
vlow vmid)))
- (list 'vec mid vmid))
-)
+ (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)
(set (nth 2 (aref math-root-vars m)) (car p)))
(setq expr-val (math-evaluate-expr expr)
jacob-val (math-evaluate-expr jacob))
- (or (and (math-constp expr-val)
- (math-constp jacob-val))
- (math-reject-arg guess "*Newton's method encountered a singularity"))
+ (unless (and (math-constp expr-val)
+ (math-constp jacob-val))
+ (math-reject-arg guess "*Newton's method encountered a singularity"))
(setq next (math-add guess (math-div (math-float (math-neg expr-val))
(math-float jacob-val)))
p guess p2 next)
limit)
(math-newton-multi expr jacob n next orig-guess limit)
(math-reject-arg nil "*Newton's method failed to converge"))
- (list 'vec next expr-val)))
-)
+ (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)
(var-DUMMY nil)
(jacob (list 'vec))
p p2 m row)
- (or (eq (car-safe var) 'vec)
- (math-reject-arg var 'vectorp))
- (or (= (length var) (1+ n))
- (math-dimension-error))
+ (unless (eq (car-safe var) 'vec)
+ (math-reject-arg var 'vectorp))
+ (unless (= (length var) (1+ n))
+ (math-dimension-error))
(setq expr (copy-sequence expr))
(while (>= n (length math-root-vars))
(let ((symb (intern (concat "math-root-v"
(while (setq p2 (cdr p2))
(setcar p2 (math-expr-subst (car p2) (car p)
(aref math-root-vars m)))))
- (or (eq (car-safe guess) 'vec)
- (math-reject-arg guess 'vectorp))
- (or (= (length guess) (1+ n))
- (math-dimension-error))
+ (unless (eq (car-safe guess) 'vec)
+ (math-reject-arg guess 'vectorp))
+ (unless (= (length guess) (1+ n))
+ (math-dimension-error))
(setq guess (copy-sequence guess)
p guess)
(while (setq p (cdr p))
(setq m (math-abs-approx guess))
(math-newton-multi expr jacob n guess guess
(if (math-zerop m) '(float 1 3) (math-mul m 10))))
- (or (eq (car-safe var) 'var)
- (math-reject-arg var "*Expected a variable"))
- (or (math-expr-contains expr var)
- (math-reject-arg expr "*Formula does not contain specified variable"))
+ (unless (eq (car-safe var) 'var)
+ (math-reject-arg var "*Expected a variable"))
+ (unless (math-expr-contains expr var)
+ (math-reject-arg expr "*Formula does not contain specified variable"))
(if (assq (car expr) calc-tweak-eqn-table)
(setq expr (math-sub (nth 1 expr) (nth 2 expr))))
(math-with-extra-prec 2
var-DUMMY guess
vlow (math-evaluate-expr expr)
vhigh vlow
- root-widen 'point)
+ math-root-widen 'point)
(if (eq (car guess) 'intv)
(progn
(or (math-constp guess) (math-reject-arg guess 'constp))
(not (Math-numberp vlow))
(not (Math-numberp vhigh)))
(math-search-root expr deriv low vlow high vhigh)
- (math-bisect-root expr low vlow high vhigh)))))))))
-)
+ (math-bisect-root expr low vlow high vhigh))))))))))
(defun calcFunc-root (expr var guess)
- (math-find-root expr var guess nil)
-)
+ (math-find-root expr var guess nil))
(defun calcFunc-wroot (expr var guess)
- (math-find-root expr var guess t)
-)
+ (math-find-root expr var guess t))
;;; 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))
(math-float a)
(if (eq (car a) 'float)
a
- (math-reject-arg a 'realp)))
-)
+ (math-reject-arg a 'realp))))
+(defvar math-min-or-max "minimum")
;;; A bracket for a minimum is a < b < c where f(b) < f(a) and f(b) < f(c).
c u vc vu))
(if (math-lessp-float a c)
(list a va b vb c vc)
- (list c vc b vb a va)))
-)
+ (list c vc b vb a va))))
(defun math-narrow-min (expr a c intv)
(let ((xvals (list a c))
(and (not yvals)
(list (nth 3 intv) min)))))
(math-reject-arg nil (format "*Unable to find a %s in the interval"
- math-min-or-max)))))
-)
+ math-min-or-max))))))
;;; "brent"
(defun math-brent-min (expr prec a va x vx b vb)
(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))
(setq v w vv vw
w x vw vx
x u vx vu)))
- (list 'vec x vx))
-)
+ (list 'vec x vx)))
;;; "powell"
(defun math-powell-min (expr n guesses prec)
(while (<= (setq i (1+ i)) n)
(setcar (nthcdr ibig (nth i xi))
(nth i (nth 1 res)))))))
- (list 'vec p fret))
-)
+ (list 'vec p fret)))
(defun math-line-min-func (expr n)
(let ((m -1))
'(var DUMMY var-DUMMY)
(list 'calcFunc-mrow '(var line-xi line-xi) (1+ m)))
(list 'calcFunc-mrow '(var line-p line-p) (1+ m)))))
- (math-evaluate-expr expr))
-)
+ (math-evaluate-expr expr)))
(defun math-line-min (f1dim line-p line-xi n prec)
(let* ((var-DUMMY nil)
(params (math-widen-min expr '(float 0 0) '(float 1 0)))
(res (apply 'math-brent-min expr prec params))
(xi (math-mul (nth 1 res) line-xi)))
- (list (math-add line-p xi) xi (nth 2 res)))
-)
+ (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)
(math-dimension-error))
(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"))
(setq guesses (cdr guesses)))
(if isvec
(list 'vec vec (nth 2 res))
- (list 'vec (nth 1 vec) (nth 2 res))))))
-)
-(setq math-min-or-max "minimum")
+ (list 'vec (nth 1 vec) (nth 2 res)))))))
(defun calcFunc-minimize (expr var guess)
(let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
(math-min-or-max "minimum"))
(math-find-minimum (math-normalize expr)
(math-normalize var)
- (math-normalize guess) nil))
-)
+ (math-normalize guess) nil)))
(defun calcFunc-wminimize (expr var guess)
(let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
(math-min-or-max "minimum"))
(math-find-minimum (math-normalize expr)
(math-normalize var)
- (math-normalize guess) t))
-)
+ (math-normalize guess) t)))
(defun calcFunc-maximize (expr var guess)
(let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
(res (math-find-minimum (math-normalize (math-neg expr))
(math-normalize var)
(math-normalize guess) nil)))
- (list 'vec (nth 1 res) (math-neg (nth 2 res))))
-)
+ (list 'vec (nth 1 res) (math-neg (nth 2 res)))))
(defun calcFunc-wmaximize (expr var guess)
(let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
(res (math-find-minimum (math-normalize (math-neg expr))
(math-normalize var)
(math-normalize guess) t)))
- (list 'vec (nth 1 res) (math-neg (nth 2 res))))
-)
+ (list 'vec (nth 1 res) (math-neg (nth 2 res)))))
(or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
(math-with-extra-prec 2
(cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
- nil))))
-)
+ nil)))))
(put 'calcFunc-polint 'math-expandable t)
(or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
(math-with-extra-prec 2
(cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
- (cdr (cdr (cdr (nth 1 data))))))))
-)
+ (cdr (cdr (cdr (nth 1 data)))))))))
(put 'calcFunc-ratint 'math-expandable t)
(setq ns (1- ns)
dy (nth ns d)))
(setq y (math-add y dy)))
- (list y dy)))
-)
+ (list y dy))))
(math-ninteg-romberg
'math-ninteg-midpoint expr
(math-float lo) (math-float hi) nil))))
- sum))
-)
+ sum)))
;;; Open Romberg method; "qromo" in section 4.4.
-(defun math-ninteg-romberg (func expr lo hi mode)
+
+;; 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)
(s nil)
(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))
(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)))
(setq curh (math-div-float curh '(float 9 0))))
ss
- (math-reject-arg nil (format "*Integral failed to converge")))))
-)
+ (math-reject-arg nil (format "*Integral failed to converge"))))))
(defun math-ninteg-evaluate (expr x mode)
(math-reject-arg res "*Integrand does not evaluate to a number"))
(if (eq mode 'inf)
(setq res (math-mul res (math-sqr x))))
- res)
-)
+ res))
-(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)))
(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)
expr (math-add x del2)
mode)))
x (math-add x del3)))
- (setq integ-temp (list it3
- (math-add (math-div (nth 1 integ-temp)
- '(float 3 0))
- (math-mul sum del)))))
- (setq integ-temp (list 1 (math-mul
- (math-sub hi lo)
- (math-ninteg-evaluate
- expr
- (math-mul (math-add lo hi) '(float 5 -1))
- mode)))))
- (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 math-ninteg-temp (list 1 (math-mul
+ (math-sub hi lo)
+ (math-ninteg-evaluate
+ expr
+ (math-mul (math-add lo hi) '(float 5 -1))
+ mode)))))
+ (nth 1 math-ninteg-temp))
;;; The following algorithms come from Numerical Recipes, chapter 14.
-(setq math-dummy-vars [(var DUMMY var-DUMMY)])
-(setq math-dummy-counter 0)
-
+(defvar math-dummy-vars [(var DUMMY var-DUMMY)])
+(defvar math-dummy-counter 0)
(defun math-dummy-variable ()
(if (= math-dummy-counter (length math-dummy-vars))
(let ((symb (intern (format "math-dummy-%d" math-dummy-counter))))
(set (nth 2 (aref math-dummy-vars math-dummy-counter)) nil)
(prog1
(aref math-dummy-vars math-dummy-counter)
- (setq math-dummy-counter (1+ math-dummy-counter)))
-)
-
+ (setq math-dummy-counter (1+ math-dummy-counter))))
+(defvar math-in-fit 0)
+(defvar calc-fit-to-trail nil)
(defun calcFunc-fit (expr vars &optional coefs data)
(let ((math-in-fit 10))
(math-with-extra-prec 2
- (math-general-fit expr vars coefs data nil)))
-)
+ (math-general-fit expr vars coefs data nil))))
(defun calcFunc-efit (expr vars &optional coefs data)
(let ((math-in-fit 10))
(math-with-extra-prec 2
- (math-general-fit expr vars coefs data 'sdev)))
-)
+ (math-general-fit expr vars coefs data 'sdev))))
(defun calcFunc-xfit (expr vars &optional coefs data)
(let ((math-in-fit 10))
(math-with-extra-prec 2
- (math-general-fit expr vars coefs data 'full)))
-)
+ (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)
- first-var
+ (math-fit-first-coef math-dummy-counter)
+ math-fit-first-var
(plain-expr expr)
orig-expr
have-sdevs need-chisq chisq
(y-filter nil)
y-dummy
(coef-filters nil)
- new-coefs
+ math-fit-new-coefs
(xy-values nil)
(weights nil)
(var-YVAL nil) (var-YVALX nil)
(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)
(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)
(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.
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)
(math-sqrt (math-mul (nth (setq j (1+ j))
(car (setq cp (cdr cp))))
sigdat))))))
- (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))
(setq vec (cons (list 'calcFunc-fitparam n) vec)
n (1- n)))
vec)
- (append (cdr new-coefs) (cdr vars))))
+ (append (cdr math-fit-new-coefs) (cdr vars))))
;; Package the result.
(math-normalize
(if (and have-sdevs (> n mm))
(list 'calcFunc-utpc chisq (- n mm))
'(var nan var-nan)))
- expr)))
-)
+ expr))))
-(setq math-in-fit 0)
-(setq calc-fit-to-trail nil)
(defun calcFunc-fitvar (x)
(if (>= math-in-fit 2)
(progn
- (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))
-)
+ (math-reject-arg x)))
(defun calcFunc-fitparam (x)
(if (>= math-in-fit 2)
(progn
- (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))
-)
+ (math-reject-arg x)))
(defun calcFunc-fitdummy (x)
(if (= math-in-fit 3)
- (nth x new-coefs)
- (math-reject-arg x))
-)
+ (nth x math-fit-new-coefs)
+ (math-reject-arg x)))
(defun calcFunc-hasfitvars (expr)
(if (Math-primp expr)
0
(if (eq (car expr) 'calcFunc-fitvar)
(nth 1 expr)
- (apply 'max (mapcar 'calcFunc-hasfitvars (cdr expr)))))
-)
+ (apply 'max (mapcar 'calcFunc-hasfitvars (cdr expr))))))
(defun calcFunc-hasfitparams (expr)
(if (Math-primp expr)
0
(if (eq (car expr) 'calcFunc-fitparam)
(nth 1 expr)
- (apply 'max (mapcar 'calcFunc-hasfitparams (cdr expr)))))
-)
+ (apply 'max (mapcar 'calcFunc-hasfitparams (cdr expr))))))
(defun math-all-vars-but (expr but)
(setq vars (delq (assoc (car-safe p) vars) vars)
p (cdr p)))
(sort (mapcar 'car vars)
- (function (lambda (x y) (string< (nth 1 x) (nth 1 y))))))
-)
+ (function (lambda (x y) (string< (nth 1 x) (nth 1 y)))))))
+
+;; The variables math-all-vars-vars (the vars for math-all-vars) and
+;; math-all-vars-found are local to math-all-vars-in, but are used by
+;; math-all-vars-rec which is called by math-all-vars-in.
+(defvar math-all-vars-vars)
+(defvar math-all-vars-found)
(defun math-all-vars-in (expr)
- (let ((vars nil)
- found)
+ (let ((math-all-vars-vars nil)
+ math-all-vars-found)
(math-all-vars-rec expr)
- vars)
-)
+ math-all-vars-vars))
(defun math-all-vars-rec (expr)
(if (Math-primp expr)
(if (eq (car-safe expr) 'var)
(or (math-const-var expr)
- (if (setq found (assoc expr vars))
- (setcdr found (1+ (cdr found)))
- (setq vars (cons (cons expr 1) vars)))))
+ (if (setq math-all-vars-found (assoc expr math-all-vars-vars))
+ (setcdr math-all-vars-found (1+ (cdr math-all-vars-found)))
+ (setq math-all-vars-vars (cons (cons expr 1) math-all-vars-vars)))))
(while (setq expr (cdr expr))
- (math-all-vars-rec (car expr))))
-)
-
-
+ (math-all-vars-rec (car expr)))))
+(provide 'calcalg3)
+;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6
+;;; calcalg3.el ends here