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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
include <mach.h>
# IN_DEVIANT -- Find deviant points with large residuals from the fit
# and reject them from the fit. The sigma of the fit residuals is calculated.
# The rejection thresholds are set at (+/-)reject*sigma. Points outside the
# rejection threshold are recorded in the reject array.
procedure in_deviantd (nl, x, y, w, rejpts, npts, nvars, low_reject,
high_reject, grow, nreject, newreject)
pointer nl # NLFIT descriptor
double x[ARB] # Input ordinates (npts * nvars)
double y[npts] # Input data values
double w[npts] # Weights
int rejpts[npts] # Points rejected
int npts # Number of input points
int nvars # Number of input variables
double low_reject, high_reject # Rejection thresholds
double grow # Rejection radius
int nreject # Number of points rejected (output)
int newreject # Number of new points rej. (output)
int i, j, i_min, i_max, ilast
double sigma, low_cut, high_cut, residual
pointer sp, residuals
begin
# # Debug.
# call eprintf (
# "in_deviant: nl=%d, npts=%d, nvars=%d, low=%g, high=%g, grow=%g\n")
# call pargi (nl)
# call pargi (npts)
# call pargi (nvars)
# call parg$t (low_reject)
# call parg$t (high_reject)
# call parg$t (grow)
# Initialize.
nreject = 0
newreject = 0
# If low_reject and high_reject are zero then just return.
if ((low_reject == double (0.0)) && (high_reject == double (0.0)))
return
# Allocate memory for the residuals.
call smark (sp)
call salloc (residuals, npts, TY_DOUBLE)
# Compute the residuals.
call nlvectord (nl, x, Memd[residuals], npts, nvars)
call asubd (y, Memd[residuals], Memd[residuals], npts)
# Compute the sigma of the residuals.
j = 0
sigma = double (0.0)
do i = 1, npts {
if ((w[i] != double (0.0)) && (rejpts[i] == NO)) {
sigma = sigma + Memd[residuals+i-1] ** 2
j = j + 1
} else if (rejpts[i] == YES)
nreject = nreject + 1
}
# If there are less than five points for the sigma calculation,
# just return.
if (j < 5) {
call sfree (sp)
return
} else
sigma = sqrt (sigma / j)
# Set the lower and upper cut limits according the the sigma value.
if (low_reject > double (0.0))
low_cut = -low_reject * sigma
else
low_cut = -MAX_REAL
if (high_reject > double (0.0))
high_cut = high_reject * sigma
else
high_cut = MAX_REAL
# Reject the residuals exceeding the rejection limits.
# A for loop is used instead of do because with region
# growing we want to modify the loop index.
for (i = 1; i <= npts; i = i + 1) {
# Do not process points with zero weigth or already rejected.
if ((w[i] == double (0.0)) || (rejpts[i] == YES))
next
# Reject point, and all other points closer than the growing
# factor.
residual = Memd[residuals + i - 1]
if ((residual > high_cut) || (residual < low_cut)) {
# Determine region to reject.
i_min = max (1, int (i - grow))
i_max = min (npts, int (i + grow))
# Reject points from the fit and flag them.
do j = i_min, i_max {
if ((abs (x[i] - x[j]) <= grow) && (w[j] != double (0.0)) &&
(rejpts[j] == NO)) {
rejpts[j] = YES
newreject = newreject + 1
ilast = j
}
}
i = ilast
}
}
# Free memory.
call sfree (sp)
end
|
