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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
include <error.h>
include "../lib/fitparams.h"
include "../lib/parser.h"
# FT_INDEF - Set zero weight for all undefined input data.
procedure ft_indef (sym, otable, rtable, wtable)
int sym # equation symbol
pointer otable # 2d observation table (modified)
pointer rtable # 1d reference table (modified)
pointer wtable # 1d weight table (modified)
#bool clgetb()
begin
# Debug ?
#if (clgetb ("debug.fitcode")) {
#call eprintf (
#"ft_indef: (sym=%d) (otable=%d) (rtable=%d) (wtable=%d)\n")
#call pargi (sym)
#call pargi (otable)
#call pargi (rtable)
#call pargi (wtable)
#}
# Check for INDEF values in reference table.
call ft_indefr (rtable, wtable)
# Check for INDEF values in fitting equation.
call ft_indefo (sym, otable, wtable)
end
# FT_INDEFR - Check reference table for INDEF values. If an INDEF value is
# found, its corresponding weight (in the weight table) is set to zero, and
# the INDEF value (in the refence table) replaced by a more suitable one.
# The latter is because the INLFIT package does not handle INDEF values at
# all, and it's better to feed it with reasonable values to avoid an overflow
# or underflow condition.
procedure ft_indefr (rtable, wtable)
pointer rtable # reference table
pointer wtable # weight table (modified)
int npts # number of points
int n
real rval, replace
#bool clgetb()
int mct_nrows()
real mct_getr()
begin
# Debug ?
#if (clgetb ("debug.fitcode")) {
#call eprintf ("ft_indefr: (rtable=%d) (wtable=%d)\n")
#call pargi (rtable)
#call pargi (wtable)
#}
# Get the number of points.
npts = mct_nrows (rtable)
# Initialize replace value to first non-INDEF value,
# if any. Otherwise set it t zero.
replace = 0.0
do n = 1, npts {
rval = mct_getr (rtable, n, 1)
if (!IS_INDEFR (rval)) {
replace = rval
break
}
}
# Loop over all data in the table.
do n = 1, npts {
# Replace values if is INDEF. Otherwise just
# update the replace value.
rval = mct_getr (rtable, n, 1)
if (IS_INDEFR (rval)) {
call mct_putr (wtable, n, 1, 0.0)
call mct_putr (rtable, n, 1, replace)
} else
replace = rval
}
# Debug ?
#call dg_dweights ("from ft_indefr", wtable)
#call dg_dref ("from ft_indefr", rtable)
end
# FT_INDEFO - Check fitting equation for INDEF values. If an INDEF value is
# found, its corresponding weight (in the weight table) is set to zero.
# Undefined values in the table are set to more suitable values, so there
# won't be problems when plotting data.
procedure ft_indefo (sym, otable, wtable)
int sym # equation symbol
pointer otable # observation table (modified)
pointer wtable # weight table (modified)
int i, n
int npts # number of points
int nvars # number of variables
real rval
pointer code # fitting equation code
pointer parval # parameter values
pointer replace # replace values
pointer sp
#bool clgetb()
int mct_nrows(), mct_maxcol()
real mct_getr()
real pr_eval()
pointer mct_getrow()
pointer pr_gsymp()
begin
# Debug ?
#if (clgetb ("debug.fitcode")) {
#call eprintf ("ft_indef: (sym=%d) (otable=%d) (wtable=%d)\n")
#call pargi (sym)
#call pargi (otable)
#call pargi (wtable)
#}
# Get the number of variables and points.
npts = mct_nrows (otable)
nvars = mct_maxcol (otable)
# Allocate space for replace values.
call smark (sp)
call salloc (replace, nvars, TY_REAL)
# Initialize replace values to first non-undefined
# value, if any. Otherwise set it t zero.
call aclrr (Memr[replace], nvars)
do i = 1, nvars {
do n = 1, npts {
rval = mct_getr (otable, n, i)
if (!IS_INDEFR (rval)) {
Memr[replace + i - 1] = rval
break
}
}
}
# Get the parameter values, and equation code.
parval = pr_gsymp (sym, PTEQSPARVAL)
code = pr_gsymp (sym, PTEQRPNFIT)
# Iterate over all the observations.
do n = 1, npts {
# Evaluate fitting equation.
rval = pr_eval (code, Memr[mct_getrow (otable, n)], Memr[parval])
# Substitute weight.
if (IS_INDEFR (rval))
call mct_putr (wtable, n, 1, 0.0)
# Substitude undefined variable values.
do i = 1, nvars {
rval = mct_getr (otable, n, i)
if (IS_INDEFR (rval))
call mct_putr (otable, n, i, Memr[replace + i - 1])
else
Memr[replace + i - 1] = rval
}
}
# Free memory.
call sfree (sp)
# Debug ?
#call dg_dweights ("from ft_indefo", wtable)
#call dg_dcatobs ("from ft_indefo", otable)
end
|
