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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
include <math/nlfit.h>
include <pkg/inlfit.h>
# IN_FIT -- Fit a function using non-linear least squares. The function
# can have an arbitrary number of independent variables. This is the main
# entry point for the non-interactive part of the INLFIT package.
procedure in_fit$t (in, nl, x, y, wts, npts, nvars, wtflag, stat)
pointer in # INLFIT pointer
pointer nl # NLFIT pointer
PIXEL x[ARB] # Ordinates (npts * nvars)
PIXEL y[npts] # Data to be fit
PIXEL wts[npts] # Weights
int npts # Number of points
int nvars # Number of variables
int wtflag # Type of weighting
int stat # Error code (output)
int i, ndeleted
pointer sp, wts1, str
int in_geti()
PIXEL in_get$t
begin
# # Debug.
# call eprintf ("in_fit: in=%d, nl=%d, npts=%d, nvars=%d\n")
# call pargi (in)
# call pargi (nl)
# call pargi (npts)
# call pargi (nvars)
# Allocate string, and rejection weight space. The latter are
# are used to mark rejected points with a zero weight before
# calling NLFIT.
call smark (sp)
call salloc (str, SZ_LINE, TY_CHAR)
call salloc (wts1, npts, TY_PIXEL)
call amov$t (wts, Mem$t[wts1], npts)
# Initialize rejected point list, and the buffer containing
# the minimum and maximum variable values.
call in_bfinit$t (in, npts, nvars)
# Set independent variable limits.
call in_limit$t (in, x, npts, nvars)
# Reinitialize.
call in_nlinit$t (in, nl)
# Check number of data points. If no points are present
# set the error flag to the appropiate value, and return.
if (npts == 0) {
stat = NO_DEG_FREEDOM
call in_puti (in, INLFITERROR, NO_DEG_FREEDOM)
call sfree (sp)
return
}
# Check the number of deleted points.
ndeleted = 0
do i = 1, npts {
if (wts[i] <= PIXEL(0.0))
ndeleted = ndeleted + 1
}
if ((npts - ndeleted) < in_geti (in, INLNFPARAMS)) {
stat = NO_DEG_FREEDOM
call in_puti (in, INLFITERROR, NO_DEG_FREEDOM)
call sfree (sp)
return
}
# Call NLFIT.
call nlfit$t (nl, x, y, wts, npts, nvars, wtflag, stat)
# Update fit status into the INLFIT structure.
call in_puti (in, INLFITERROR, stat)
# Do pixel rejection and refit, if at least one of the rejection
# limits is positive. Otherwise clear number of rejected points.
if (in_get$t (in, INLLOW) > PIXEL (0.0) ||
in_get$t (in, INLHIGH) > PIXEL (0.0)) {
call in_reject$t (in, nl, x, y, Mem$t[wts1], npts, nvars, wtflag)
if (in_geti (in, INLNREJPTS) > 0) {
do i = 1, npts {
if (Mem$t[wts1+i-1] > PIXEL(0.0))
wts[i] = Mem$t[wts1+i-1]
}
}
stat = in_geti (in, INLFITERROR)
} else
call in_puti (in, INLNREJPTS, 0)
# Free memory.
call sfree (sp)
end
|
