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
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
|
include "../lib/parser.h"
define DET_TOL 1.0e-20
# PH_IEREQN -- Compute error estimates for the photometric indices based
# on the value of any user-supplied errors equations. Use the chi-fitting
# method developed by Bevington.
int procedure ph_iereqn (params, a, nobs, deltaa, aindex, errors, nstd,
seteqn, setvals, serrors, nset, nustd, eqindex, neq)
pointer params[ARB] # input array of pointers to the fitted parameters
real a[ARB] # array of observed and catalog variables
int nobs # the number of observed variables
real deltaa[ARB] # array of increments for the catalog variables
int aindex[ARB] # list of active catalog variables
real errors[ARB] # output array of error estimates
int nstd # number of catalog variables to be fit
int seteqn[ARB] # array of set equation codes
real setvals[ARB] # the set equation values
real serrors[ARB] # output array of set equation errors
int nset # the number of set equations
int nustd # the number of active catalog variables
int eqindex[ARB] # array of equations indices
int neq # number of equations
int i, j, sym, nerrors
pointer errcode
real rms, det
int pr_gsym()
pointer pr_gsymp()
real pr_eval(), ph_accum(), ph_incrms()
include "invert.com"
begin
# Evaluate the reference and error equations.
nerrors = 0
do i = 1, neq {
sym = pr_gsym (eqindex[i], PTY_TRNEQ)
Memr[py+i-1] = pr_eval (pr_gsymp (sym, PTEQRPNREF), a,
Memr[params[eqindex[i]]])
if (IS_INDEFR(Memr[py+i-1]))
return (0)
errcode = pr_gsymp (sym, PTEQRPNERROR)
if (errcode == NULL)
Memr[pyerr+i-1] = 0.0
else {
Memr[pyerr+i-1] = pr_eval (errcode, a,
Memr[params[eqindex[i]]])
if (IS_INDEFR(Memr[pyerr+i-1]))
Memr[pyerr+i-1] = 0.0
else
nerrors = nerrors + 1
}
}
# Return if there are no error equations.
if (nerrors <= 0)
return (0)
# Compute each equations contribution to the to the total error.
call aclrr (errors, nstd)
call aclrr (serrors, nset)
do i = 1, neq {
# Add in the appropriate error term.
if (Memr[pyerr+i-1] <= 0.0)
next
call amovr (a[nobs+1], Memr[pafit], nstd)
Memr[py+i-1] = Memr[py+i-1] + Memr[pyerr+i-1]
# Accumulate the matrices and vectors, do the inversion, and
# compute the new parameter increments.
rms = ph_accum (params, Memr[py], Memr[pyfit], eqindex, neq, a,
nobs, deltaa, aindex, Memr[pda], nstd, Memd[palpha],
Memr[pbeta], Memi[pik], Memi[pjk], nustd, det)
# Compute the contribution to the sum of the squares of the errors.
#if ((! IS_INDEFR(rms)) && (abs (det) >= DET_TOL)) {
if (! IS_INDEFR(rms)) {
rms = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, neq,
a, nobs, Memr[pda], aindex, nstd, rms)
do j = 1, nstd {
if (aindex[j] == 0)
next
errors[j] = errors[j] + (a[nobs+j] - Memr[pafit+j-1]) ** 2
}
do j = 1, nset {
serrors[j] = serrors[j] + (pr_eval (seteqn[j], a,
Memr[params[1]]) - setvals[j]) ** 2
}
}
# Reset the error term.
Memr[py+i-1] = Memr[py+i-1] - Memr[pyerr+i-1]
call amovr (Memr[pafit], a[nobs+1], nstd)
}
# Compute the errors themselves.
do i = 1, nstd {
if (errors[i] <= 0.0)
errors[i] = 0.0
else
errors[i] = sqrt (errors[i])
}
do i = 1, nset {
if (serrors[i] <= 0.0)
serrors[i] = 0.0
else
serrors[i] = sqrt (serrors[i])
}
return (nerrors)
end
# PH_IERVAL -- Compute error estimates for the photometric indices based
# on the value of any user-supplied errors equations.
int procedure ph_ierval (params, a, eindex, nobs, deltaa, aindex, errors,
nstd, seteqn, setvals, serrors, nset, nustd, eqindex, neq)
pointer params[ARB] # input array of pointers to the fitted parameters
real a[ARB] # array of observed and catalog variables
int eindex[ARB] # indices of observed and catalog variables with errors
int nobs # the number of observed variables
real deltaa[ARB] # array of increments for the catalog variables
int aindex[ARB] # list of active catalog variables
real errors[ARB] # output array of error estimates
int nstd # number of catalog variables to be fit
int seteqn[ARB] # array of set equation codes
real setvals[ARB] # the set equation values
real serrors[ARB] # output array of set equation errors
int nset # the number of set equations
int nustd # the number of active catalog variables
int eqindex[ARB] # array of equation indices
int neq # number of equations
int i, j, k, sym, nerrors
real atemp, rms, det
int pr_gsym()
pointer pr_gsymp()
real pr_eval(), ph_accum(), ph_incrms()
include "invert.com"
begin
# Initialize.
nerrors = 0
call aclrr (errors, nstd)
call aclrr (serrors, nset)
# Loop over the observational variables.
do j = 1, nobs {
# Use only variables with errors.
if (eindex[j] <= 0)
next
nerrors = nerrors + 1
if (IS_INDEFR(a[eindex[j]]) || (a[eindex[j]] <= 0.0))
next
atemp = a[j]
a[j] = a[j] + a[eindex[j]]
call amovr (a[nobs+1], Memr[pafit], nstd)
# Evaluate the reference equations.
do i = 1, neq {
sym = pr_gsym (eqindex[i], PTY_TRNEQ)
Memr[py+i-1] = pr_eval (pr_gsymp (sym, PTEQRPNREF), a,
Memr[params[eqindex[i]]])
if (IS_INDEFR(Memr[py+i-1]))
return (0)
}
# Accumulate the matrices and vectors, do the inversion, and
# compute the new parameter increments.
#call eprintf ("errors before accum\n")
rms = ph_accum (params, Memr[py], Memr[pyfit], eqindex, neq, a,
nobs, deltaa, aindex, Memr[pda], nstd, Memd[palpha],
Memr[pbeta], Memi[pik], Memi[pjk], nustd, det)
#call eprintf ("errors after accum\n")
# Compute the contribution to the sum of the squares of the errors.
#if ((! IS_INDEFR(rms)) && (abs (det) >= DET_TOL)) {
if (! IS_INDEFR(rms)) {
#call eprintf ("errors before accum\n")
rms = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, neq,
a, nobs, Memr[pda], aindex, nstd, rms)
#call eprintf ("errors before accum\n")
do k = 1, nstd {
if (aindex[k] <= 0)
next
errors[k] = errors[k] + (a[nobs+k] - Memr[pafit+k-1]) ** 2
}
a[j] = atemp
do k = 1, nset {
serrors[k] = serrors[k] + (pr_eval (seteqn[k], a,
Memr[params[1]]) - setvals[k]) ** 2
}
call amovr (Memr[pafit], a[nobs+1], nstd)
} else {
a[j] = atemp
call amovr (Memr[pafit], a[nobs+1], nstd)
}
}
# Compute the errors themselves.
do i = 1, nstd {
if (errors[i] <= 0.0)
errors[i] = 0.0
else
errors[i] = sqrt (errors[i])
}
do i = 1, nset {
if (serrors[i] <= 0.0)
serrors[i] = 0.0
else
serrors[i] = sqrt (serrors[i])
}
return (nerrors)
end
# PH_ERVAL -- Compute the error in an equation by summing the effects of errors
# in the observational variables.
real procedure ph_erval (symfit, fitval, params, a, aindex, eindex, nobsvars)
pointer symfit # pointer to the fit equation
real fitval # the best fit value
real params[ARB] # the fitted equation parameters
real a[ARB] # the observed and catalog variable values
int aindex[ARB] # the list of active variables
int eindex[ARB] # the list
int nobsvars # the number of observed variables
int i, nerrors
real errval, afit
real pr_eval()
begin
# Initialize.
nerrors = 0
errval = 0.0
# Loop over the observed variables.
do i = 1, nobsvars {
# Skip observed variables with no errors.
if (eindex[i] <= 0)
next
nerrors = nerrors + 1
if ((IS_INDEFR(a[eindex[i]])) || (a[eindex[i]] <= 0.0))
next
# Evaluate the contribution of each observed variable to the
# total error.
afit = a[i]
a[i] = a[i] + a[eindex[i]]
errval = errval + (pr_eval (symfit, a, params) - fitval) ** 2
a[i] = afit
}
if (nerrors <= 0)
return (INDEFR)
else
return (sqrt (errval))
end
|
