1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
|
include <mach.h>
include <math/nlfit.h>
$if (datatype == r)
include "nlfitdefr.h"
$else
include "nlfitdefd.h"
$endif
# NLITER - Routine to compute one iteration of the fitting process
procedure nliter$t (nl, x, z, w, npts, nvars, ier)
pointer nl # pointer to nl fitting structure
PIXEL x[ARB] # independent variables (npts * nvars)
PIXEL z[ARB] # function values (npts)
PIXEL w[ARB] # weights (npts)
int npts # number of points
int nvars # number of independent variables
int ier # error code
int i, index
PIXEL nlacpts$t(), nlresid$t()
begin
# Do the initial accumulation.
NL_OLDSQ(nl) = nlacpts$t (nl, x, z, w, npts, nvars)
# Set up temporary parameter array.
call amov$t (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
repeat {
# Solve the matrix.
call nlsolve$t (nl, ier)
if (ier == NO_DEG_FREEDOM)
return
# Increment the parameters.
do i = 1, NL_NFPARAMS(nl) {
index = PLIST(NL_PLIST(nl)+i-1)
TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
DPARAM(NL_DPARAM(nl)+i-1)
}
# Reaccumulate the residuals and increment lambda.
NL_SUMSQ(nl) = nlresid$t (nl, x, z, w, npts, nvars)
#if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
call amov$t (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
NL_NPARAMS(nl))
NL_LAMBDA(nl) = PIXEL (0.10) * NL_LAMBDA(nl)
break
} else
NL_LAMBDA(nl) = min (PIXEL(LAMBDAMAX),
PIXEL (10.0) * NL_LAMBDA(nl))
} until (NL_LAMBDA(nl) <= PIXEL(0.0) ||
NL_LAMBDA(nl) >= PIXEL(LAMBDAMAX))
end
|
