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