diff options
Diffstat (limited to 'noao/digiphot/photcal/evaluate')
| -rw-r--r-- | noao/digiphot/photcal/evaluate/README | 3 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/invert.com | 15 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/mkpkg | 18 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/phcheck.x | 240 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/pherrors.x | 278 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/phinvert.x | 619 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/phminv.f | 114 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/phprint.x | 110 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/t_evalfit.x | 495 | ||||
| -rw-r--r-- | noao/digiphot/photcal/evaluate/t_invertfit.x | 792 |
10 files changed, 2684 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/evaluate/README b/noao/digiphot/photcal/evaluate/README new file mode 100644 index 00000000..879dfcee --- /dev/null +++ b/noao/digiphot/photcal/evaluate/README @@ -0,0 +1,3 @@ +This subdirectory contains the code for evaluatin the fit found by FITCOEFFS +either in the forward sense using the EVALUATE task or in the inverse sense +using the INVEVAL task. diff --git a/noao/digiphot/photcal/evaluate/invert.com b/noao/digiphot/photcal/evaluate/invert.com new file mode 100644 index 00000000..96af60d9 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/invert.com @@ -0,0 +1,15 @@ +pointer py # pointer to the reference equation values +pointer pyfit # pointer to the fit equation values +pointer pa # pointer to the parameters to be fit +pointer pdelta # pointer to the parameter increments +pointer pda # pointer to the temporary parameter increments +pointer palpha # pointer to the accumulated matrix +pointer pbeta # pointer to the accumulated right side vector +pointer pik # pointer to working index array for matrix inversion +pointer pjk # pointer to working index array for matrix inversion + +pointer pyerr # pointer to the error equation values +pointer pafit # pointer to the temporary fit array + +common /invertcom / py, pyfit, pa, pdelta, pda, palpha, pbeta, pik, pjk +common /invertcom / pyerr, pafit diff --git a/noao/digiphot/photcal/evaluate/mkpkg b/noao/digiphot/photcal/evaluate/mkpkg new file mode 100644 index 00000000..667d92e3 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/mkpkg @@ -0,0 +1,18 @@ +# The MKPKG file for the INVEVALUATE and PHOTPROC task. + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + t_invertfit.x "../lib/io.h" "../lib/parser.h" <math/nlfit.h> \ + <error.h> "../lib/preval.h" + t_evalfit.x "../lib/io.h" "../lib/parser.h" <math/nlfit.h> \ + <error.h> + phinvert.x <mach.h> "invert.com" "../lib/parser.h" + pherrors.x "../lib/parser.h" "invert.com" + phprint.x "../lib/parser.h" + phcheck.x "../lib/parser.h" "../lib/preval.h" + phminv.f + ; diff --git a/noao/digiphot/photcal/evaluate/phcheck.x b/noao/digiphot/photcal/evaluate/phcheck.x new file mode 100644 index 00000000..2f729930 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/phcheck.x @@ -0,0 +1,240 @@ +include "../lib/parser.h" +include "../lib/preval.h" + +# PH_SETEQN -- Determine whether a PHOTCAL expression includes any set +# equations which contain references to catalog variables. Return a YES +# or NO status. + +int procedure ph_seteqn (code) + +pointer code # the photcal equation code + +int ip, ins +pointer setcode +int pr_geti(), pr_gsym(), ph_cateqn() +pointer pr_gsymp() + +begin + # Return NO if there are no defined set equations. + if (pr_geti (NSETEQS) <= 0) + return (NO) + + # Initialize. + ip = 0 + ins = Memi[code+ip] + + # Loop over the equation parsing instructions. + while (ins != PEV_EOC) { + + switch (ins) { + case PEV_SETEQ: + ip = ip + 1 + setcode = pr_gsymp (pr_gsym (Memi[code+ip], PTY_SETEQ), + PSEQRPNEQ) + if (ph_cateqn (setcode) == YES) + return (YES) + case PEV_NUMBER, PEV_CATVAR, PEV_OBSVAR, PEV_PARAM: + ip = ip + 1 + case PEV_EXTEQ, PEV_TRNEQ: + ip = ip + 1 + default: + # do nothing + } + + ip = ip + 1 + ins = Memi[code+ip] + } + + return (NO) +end + + +# PH_CATEQN -- Given a PHOTCAL set equation determine whether there are any +# references in it to catalog variables. Return a YES or NO status. + +int procedure ph_cateqn (code) + +pointer code # pointer to the equation code + +int ip, ins +int pr_geti() + +begin + # Return NO if there are no catalog variables. + if (pr_geti (NCATVARS) <= 0) + return (NO) + + # Initialize. + ip = 0 + ins = Memi[code+ip] + + # Loop over the equation parsing instructions. + while (ins != PEV_EOC) { + + switch (ins) { + case PEV_CATVAR: + return (YES) + case PEV_NUMBER, PEV_OBSVAR, PEV_PARAM: + ip = ip + 1 + case PEV_SETEQ, PEV_EXTEQ, PEV_TRNEQ: + ip = ip + 1 + default: + # do nothing + } + + ip = ip + 1 + ins = Memi[code+ip] + } + + return (NO) +end + + +# PH_SETVAR -- Check to see whether the fit expression portion of a PHOTCAL +# transformation equation contains references to set equations which include +# previously unreferenced catalog variables. Return the number of previously +# unreferenced set equations. + +int procedure ph_setvar (eqn, sym, cmap, omap, tuservars, uservars, usererrs, + nstdvars, nobsvars, eqset, maxnset, userset, nset) + +int eqn # the equation number +pointer sym # pointer to the equation symbol +pointer cmap # pointer to catalog variable column map +pointer omap # pointer to observations variable column map +int tuservars[nstdvars,ARB] # list of active catalog variables per equation +int uservars[ARB] # list of active catalog variables +int usererrs[ARB] # list of active observational error columns +int nstdvars # total number of catalog variables +int nobsvars # total number of observations variables +int eqset[maxnset,ARB] # set equation table +int maxnset # maximum number of set equations +int userset[ARB] # the list of active set equations +int nset # number of set equations + +bool new +int i, ip, ins, nvar +pointer fitcode, setsym +int pr_gsym(), ph_catvar() +pointer pr_gsymp() + +begin + # Get the fit expression symbols. + fitcode = pr_gsymp (sym, PTEQRPNFIT) + + # Initialize. + nvar = 0 + ip = 0 + ins = Memi[fitcode+ip] + + # Loop over the fit expression. + while (ins != PEV_EOC) { + + switch (ins) { + case PEV_SETEQ: + ip = ip + 1 + setsym = pr_gsym (Memi[fitcode+ip], PTY_SETEQ) + if (ph_catvar (eqn, setsym, cmap, omap, tuservars, uservars, + usererrs, nstdvars, nobsvars) == YES) { + new = true + do i = 1, nset + nvar { + if (setsym != userset[i]) + next + eqset[i,eqn] = setsym + new = false + break + } + if (new) { + nvar = nvar + 1 + eqset[nset+nvar,eqn] = setsym + userset[nset+nvar] = setsym + } + } + case PEV_NUMBER, PEV_CATVAR, PEV_OBSVAR, PEV_PARAM: + ip = ip + 1 + case PEV_EXTEQ, PEV_TRNEQ: + ip = ip + 1 + default: + # do nothing + } + + ip = ip + 1 + ins = Memi[fitcode+ip] + } + + return (nvar) +end + + +# PH_CATVAR -- Determine whether any catalog variables in the set +# equation have not been previously referenced in the transformation +# equations. + +int procedure ph_catvar (eqn, setsym, cmap, omap, tuservars, uservars, + usererrs, nstdvars, nobsvars) + +int eqn # the transformation equation number +int setsym # the set equation symbol +pointer cmap # pointer to catalog variable column map +pointer omap # pointer to observations variable column map +int tuservars[nstdvars,ARB] # active variables as a function of eqn +int uservars[ARB] # currently active catalog variables +int usererrs[ARB] # currently active observations varaiables errs +int nstdvars # the total number of catalog variables +int nobsvars # the total number of observations variables + +int setcode, ip, ins, stat, vcol, ecol +pointer catsym, obsym +int pr_gsym(), pr_gsymi(), pr_findmap1() +pointer pr_gsymp() + +begin + # Get the set equation code. + setcode = pr_gsymp (setsym, PSEQRPNEQ) + + # Initialize. + ip = 0 + ins = Memi[setcode+ip] + stat = NO + + while (ins != PEV_EOC) { + + switch (ins) { + case PEV_CATVAR: + ip = ip + 1 + catsym = pr_gsym (Memi[setcode+ip] - nobsvars, PTY_CATVAR) + vcol = pr_gsymi (catsym, PINPCOL) + vcol = pr_findmap1 (cmap, vcol) + if (uservars[vcol] <= 0) { + tuservars[vcol,eqn] = -(vcol + nobsvars) + uservars[vcol] = -(vcol + nobsvars) + stat = YES + } + + case PEV_OBSVAR: + ip = ip + 1 + obsym = pr_gsym (Memi[setcode+ip], PTY_OBSVAR) + vcol = pr_gsymi (obsym, PINPCOL) + ecol = pr_gsymi (obsym, PINPERRCOL) + vcol = pr_findmap1 (omap, vcol) + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (omap, ecol) + if (! IS_INDEFI(ecol)) + usererrs[vcol] = ecol + + case PEV_NUMBER, PEV_PARAM: + ip = ip + 1 + + case PEV_SETEQ, PEV_EXTEQ, PEV_TRNEQ: + ip = ip + 1 + + default: + # do nothing + } + + ip = ip + 1 + ins = Memi[setcode+ip] + } + + return (stat) +end diff --git a/noao/digiphot/photcal/evaluate/pherrors.x b/noao/digiphot/photcal/evaluate/pherrors.x new file mode 100644 index 00000000..77a95391 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/pherrors.x @@ -0,0 +1,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 diff --git a/noao/digiphot/photcal/evaluate/phinvert.x b/noao/digiphot/photcal/evaluate/phinvert.x new file mode 100644 index 00000000..9d7ee14b --- /dev/null +++ b/noao/digiphot/photcal/evaluate/phinvert.x @@ -0,0 +1,619 @@ +include <mach.h> +include "../lib/parser.h" + +# PH_IVINIT -- Preallocate the space required by the inversion code. + +procedure ph_ivinit (nstd, nustd, neq) + +int nstd # number of catalog variables +int nustd # number of catalog variables to be fit +int neq # number of equations + +include "invert.com" + +begin + call malloc (py, neq, TY_REAL) + call malloc (pyfit, neq, TY_REAL) + call malloc (pa, nstd, TY_REAL) + call malloc (pdelta, nstd, TY_REAL) + call malloc (pda, nstd, TY_REAL) + + call malloc (palpha, nustd * nustd, TY_DOUBLE) + call malloc (pbeta, nustd, TY_REAL) + call malloc (pik, nustd, TY_INT) + call malloc (pjk, nustd, TY_INT) + + call malloc (pyerr, neq, TY_REAL) + call malloc (pafit, nstd, TY_REAL) +end + + +# PH_IVFREE -- Free the space required by the inversion code. + +procedure ph_ivfree () + +include "invert.com" + +begin + call mfree (py, TY_REAL) + call mfree (pyfit, TY_REAL) + call mfree (pa, TY_REAL) + call mfree (pdelta, TY_REAL) + call mfree (pda, TY_REAL) + + call mfree (palpha, TY_DOUBLE) + call mfree (pbeta, TY_REAL) + call mfree (pik, TY_INT) + call mfree (pjk, TY_INT) + + call mfree (pyerr, TY_REAL) + call mfree (pafit, TY_REAL) +end + + +# PH_OBJCHECK -- Check that the equations for this particular star are +# invertable. + +int procedure ph_objcheck (params, a, vartable, nstdvars, nreq, eqset, + maxnset, vindex, nvar, eqindex, neq) + +pointer params[ARB] # array of pointers to the fitted parameters +real a[ARB] # array of observed and catalog variables +int vartable[nstdvars,ARB] # table of variables as a function of equation +int nstdvars # the total number of catalog variables +int nreq # the total number of reference equations +int eqset[maxnset,ARB] # set equation table +int maxnset # maximum number of set equations +int vindex[ARB] # output index of variables +int nvar # number of variables used in the equations +int eqindex[ARB] # output index of equations +int neq # number of equations to be reduced + +int i, j, sym, ncat, nset +real rval +int pr_gsym() +pointer pr_gsymp() +real pr_eval() + +begin + # Initialize + call aclri (vindex, nstdvars) + call aclri (eqindex, nreq) + + # Evalute the reference equations. + neq = 0 + do i = 1, nreq { + sym = pr_gsym (i, PTY_TRNEQ) + rval = pr_eval (pr_gsymp (sym, PTEQRPNREF), a, Memr[params[i]]) + if (IS_INDEFR(rval)) + next + neq = neq + 1 + eqindex[neq] = i + } + + # If there is no data return. + if (neq <= 0) + return (ERR) + + # Determine which variables are used by the reduced set of equations. + do i = 1, neq { + do j = 1, nstdvars { + if (vartable[j,eqindex[i]] == 0) + next + vindex[j] = vartable[j,eqindex[i]] + } + } + + # Deterine which set equations are used by the reduced set of equations + nset = 0 + do j = 1, maxnset { + do i = 1, neq { + if (eqset[j,eqindex[i]] == 0) + next + nset = nset + 1 + break + } + } + + # Count the number of variables. + nvar = 0 + ncat = 0 + do j = 1, nstdvars { + if (vindex[j] == 0) + next + nvar = nvar + 1 + if (vindex[j] > 0) + ncat = ncat + 1 + } + + if ((ncat + nset) > neq) + return (ERR) + + return (OK) +end + + +define MAX_NITER1 10 # maximum number of iterations +define MAX_NITER2 50 # maximum number of trials +define DET_TOL 1.0E-20 # minimum value of the determinant +#define DET_TOL 0.0 # minimum value of the determinant +define MIN_DELTA 0.01 # minimum absolute value of + # parameter increments + +# PH_INVERT -- Invert the transformation to compute the standard indices. + +int procedure ph_invert (params, a, nobs, deltaa, aindex, nstd, + nustd, eqindex, nueq) + +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] # index of active catalog variables +int nstd # total number of catalog variables +int nustd # number of catalog variables to be fit +int eqindex[ARB] # the equation index +int nueq # total number of equations used + +int i, sym, niter +real stdev1, stdev2, det, rms +int pr_gsym() +pointer pr_gsymp() +real pr_eval(), ph_accum(), ph_incrms() + +include "invert.com" + +begin + # Evalute the reference equations. + do i = 1, nueq { + 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 (ERR) + } + + # Initialize the parameter increments. This will be incremented + # each time through the fitting loop. + call amovr (deltaa, Memr[pdelta], nstd) + + # Accumulate the matrices and vectors, do the inversion and compute + # the new parameter increments. The fit will terminate when the + # determinant of the inversion matrix becomes < 1.0e-20, if the + # standard deviation of the fit begins to increase, if the rms of + # the fit < EPSILONR, or the maximum number of iterations is + # exceeded, whichever comes first. + + niter = 0 + + repeat { + + # Compute the curvature currection. Return INDEFR if there are + # bad data in the fit. Terminate the fit if the determinant of + # the curvature matrix is too close to zero. + + stdev1 = ph_accum (params, Memr[py], Memr[pyfit], eqindex, nueq, + a, nobs, Memr[pdelta], aindex, Memr[pda], nstd, + Memd[palpha], Memr[pbeta], Memi[pik], Memi[pjk], nustd, det) + + #call eprintf ("acc: niter=%d det=%g stdev1=%g\n") + #call pargi (niter+1) + #call pargr (det) + #call pargr (stdev1) + + # Return if there is INDEF data in the fit. + if (IS_INDEFR(stdev1)) + return (ERR) + + # Check the size of the determinant but force at least one fit. + if ((abs (det) < DET_TOL) && (niter > 0)) + break + + # Find the new parameter values. + call amovr (a, Memr[pa], nstd) + stdev2 = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, nueq, + a, nobs, Memr[pda], aindex, nstd, stdev1) + + #call eprintf ("inc: niter=%d det=%g stdev2=%g\n") + #call pargi (niter+1) + #call pargr (det) + #call pargr (stdev2) + + # Check the new values. + if (IS_INDEFR(stdev2) || (stdev2 >= stdev1)) { + if (niter == 0) + return (ERR) + else { + call amovr (Memr[pa], a, nstd) + return (OK) + } + } + + # User the new deltas and increment the fit counters. + call anegr (Memr[pda], Memr[pdelta], nstd) + call ph_deltamin (Memr[pdelta], MIN_DELTA, Memr[pdelta], nstd) + niter = niter + 1 + rms = sqrt (stdev2) + + } until ((niter == MAX_NITER1) || (rms <= EPSILONR)) + + return (OK) +end + + +# PH_ACCUM -- Accumulate the matrix of second derivatives and vector +# of first derivatives required for parabolic expansion of the rms +# non-linear least squares fitting technique. This code is a modification +# of the Bevington CHIFIT subroutine where the reduced chi-squared has +# been replaced by the rms. The original CHIFIT cannot be used to fit the +# case when there are n data points and n unknowns since the number of +# degrees of freedom is zero and hence so is the reduced chi-squared. + +real procedure ph_accum (params, y, yfit, eqindex, neq, a, nobs, deltaa, + aindex, da, nstd, alpha, beta, ik, jk, nustd, det) + +pointer params[ARB] # input array of ptrs to the fitted parameters +real y[ARB] # array of reference equation values +real yfit[ARB] # array of fitted reference equation values +int eqindex[ARB] # the equation indices +int neq # number of equations to be inverted +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] # index of active catalog variables +real da[ARB] # new array of parameter increments +int nstd # number of catalog variables +double alpha[nustd,ARB] # the matrix of second derivatives +real beta[ARB] # the vector of first derivatives +int ik[ARB] # working array for the matrix inversion +int jk[ARB] # working array for the matrix inversion +int nustd # The number of catalog variables to be fit +real det # determinant of inverted matrix + +int i, j, k, nj, nk, sym +real aj, ak, rms1, rms2, rms3 +int pr_gsym() +pointer pr_gsymp() +real pr_eval() + +begin + # Compute the initial rms, checking for INDEF valued variables. + rms1 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + if (IS_INDEFR(yfit[i])) + return (INDEFR) + rms1 = rms1 + (y[i] - yfit[i]) ** 2 + } + rms1 = rms1 / neq + + nj = 0 + do j = 1, nstd { + + # Check the status of the parameter. + if (aindex[j] == 0) + next + nj = nj + 1 + + # Increment each parameter. + aj = a[j+nobs] + a[j+nobs] = aj + deltaa[j] + + # Compute a new rms. + rms2 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms2 = rms2 + (y[i] - yfit[i]) ** 2 + } + rms2 = rms2 / neq + + # Begin accumulating the diagonal elements . + alpha[nj,nj] = rms2 - 2.0 * rms1 + beta[nj] = -rms2 + + # Accumulate the non-diagonal elements. + nk = 0 + do k = 1, nstd { + + if (aindex[k] == 0) + next + nk = nk + 1 + + if ((nk - nj) == 0) + next + else if ((nk - nj) < 0) { + alpha[nk,nj] = (alpha[nk,nj] - rms2) / 2.0 + alpha[nj,nk] = alpha[nk,nj] + next + } + + alpha[nj,nk] = rms1 - rms2 + ak = a[k+nobs] + a[k+nobs] = ak + deltaa[k] + + rms3 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms3 = rms3 + (y[i] - yfit[i]) ** 2 + } + rms3 = rms3 / neq + + alpha[nj,nk] = alpha[nj,nk] + rms3 + a[k+nobs] = ak + } + + # Continue accumulating the diagonal elements. + a[j+nobs] = aj - deltaa[j] + rms3 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms3 = rms3 + (y[i] - yfit[i]) ** 2 + } + rms3 = rms3 / neq + + a[j+nobs] = aj + alpha[nj,nj] = (alpha[nj,nj] + rms3) / 2.0 + beta[nj] = (beta[nj] + rms3) / 4.0 + } + + # Eliminate any curvature from the matrix of second derivatives. + do j = 1, nj { + if (alpha[j,j] > 0.0) + next + if (alpha[j,j] < 0.0) + alpha[j,j] = -alpha[j,j] + else + alpha[j,j] = 0.01 + do k = 1, nk { + if ((k - j) == 0) + next + alpha[j,k] = 0.0 + alpha[k,j] = 0.0 + } + } + + # Invert the matrix. + call phminv (alpha, ik, jk, nj, det) + + # Increment the parameters. + nj = 0 + do j = 1, nstd { + da[j] = 0.0 + if (aindex[j] == 0) + next + nj = nj + 1 + nk = 0 + do k = 1, nstd { + if (aindex[k] == 0) + next + nk = nk + 1 + da[j] = da[j] + beta[nk] * alpha[nj,nk] + } + da[j] = 0.2 * da[j] * deltaa[j] + } + + # If the determinate is too small increment the parameters by + # deltas and try again. + #if (abs (det) < DET_TOL) { + #call eprintf ("using approx\n") + #do j = 1, nstd { + #if (aindex[j] == 0) { + #da[j] = 0.0 + #next + #} + #call eprintf ("i=%d dain=%g ") + #call pargi (j) + #call pargr (da[j]) + #if (da[j] > 0.0) + #da[j] = abs (deltaa[j]) + #else if (da[j] < 0.0) + #da[j] = -abs (deltaa[j]) + #else + #da[j] = 0.0 + #call eprintf ("daout=%g\n") + #call pargr (da[j]) + #} + #} + + return (rms1) +end + + +# PH_INCRMS -- Increment the parameters until three values of the rms are found +# which bracket the best fitting data point and fit a parabola to the three +# different rms points. + +real procedure ph_incrms (params, y, yfit, eqindex, neq, a, nobs, da, aindex, + nstd, rms1) + +pointer params[ARB] # input array of ptrs to the fitted parameters +real y[ARB] # the array of reference equation values +real yfit[ARB] # the array of fitted reference equation values +int eqindex[ARB] # list of active equations +int neq # the number of equations +real a[ARB] # array of observed and fitted variables +int nobs # number of observed variables +real da[ARB] # the parameter increments +int aindex[ARB] # the index of active catalog variables +int nstd # number of catalog variables +real rms1 # the first rms point + +int j, i, sym, niter +real orms2, rms2, orms3, rms3, delta, rms +int pr_gsym() +pointer pr_gsymp() +real pr_eval() + +begin + # Adjust the parameters. + rms = rms1 + do j = 1, nstd { + if (aindex[j] == 0) + next + a[j+nobs] = a[j+nobs] + da[j] + } + + # Alter the parameter increments until the rms starts to decrease. + orms2 = MAX_REAL + niter = 0 + repeat { + + rms2 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms2 = rms2 + (y[i] - yfit[i]) ** 2 + } + rms2 = rms2 / neq + #call eprintf (" niter=%d rms1=%g rms2=%g\n") + #call pargi (niter) + #call pargr (rms1) + #call pargr (rms2) + if (rms2 <= 0.0) + break + if ((rms1 - rms2) >= 0.0) + break + + # If rms2 does not decrease and does not change from one iteration + # to the next we are probably near the precision limits of the + # computer or the computed curvature has the wrong sign. In that + # case quit. + + if (orms2 < MAX_REAL) { + if (abs ((rms2 - orms2) / orms2) < EPSILONR) + return (rms) + } + + orms2 = rms2 + niter = niter + 1 + if (niter >= MAX_NITER2) + return (rms) + + do j = 1, nstd { + if (aindex[j] == 0) + next + #da[j] = da[j] / 2.0 + a[j+nobs] = a[j+nobs] - da[j] + } + + + } + + # Alter the parameter increments until the rms starts to increase. + orms3 = MAX_REAL + niter = 0 + repeat { + + do j = 1, nstd { + if (aindex[j] == 0) + next + a[j+nobs] = a[j+nobs] + da[j] + } + + rms3 = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms3 = rms3 + (y[i] - yfit[i]) ** 2 + } + rms3 = rms3 / neq + #call eprintf (" niter=%d rms1=%g rms2=%g rms3=%g\n") + #call pargi (niter) + #call pargr (rms1) + #call pargr (rms2) + #call pargr (rms3) + if ((rms3 - rms2) >= 0.0) + break + + if (orms3 < MAX_REAL) { + if (abs ((rms3 - orms3) / orms3) < EPSILONR) + return (rms) + } + + niter = niter + 1 + if (niter >= MAX_NITER2) + return (rms) + + + orms3 = rms3 + rms1 = rms2 + rms2 = rms3 + + } + + # Fit a parabola to the three values of the rms that bracket the fit. + if (rms3 <= rms2) + delta = 1.0 + else + delta = 1.0 / (1.0 + (rms1 - rms2) / (rms3 - rms2)) + 0.5 + do j = 1, nstd { + if (aindex[j] == 0) + next + a[nobs+j] = a[nobs+j] - delta * da[j] + } + + rms = 0.0 + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + rms = rms + (y[i] - yfit[i]) ** 2 + } + rms = rms / neq + + if ((rms2 - rms) < 0.0) { + do j = 1, nstd { + if (aindex[j] == 0) + next + a[nobs+j] = a[nobs+j] + (delta - 1.) * da[j] + } + do i = 1, neq { + sym = pr_gsym (eqindex[i], PTY_TRNEQ) + yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a, + Memr[params[eqindex[i]]]) + } + rms = rms2 + } + + #call eprintf (" incr: rms1=%g rms2=%g rms3=%g rms=%g\n") + #call pargr (rms1) + #call pargr (rms2) + #call pargr (rms3) + #call pargr (rms) + + return (rms) +end + + +# PH_DELTAMIN -- Check to make sure that the absolute value of the deltaa +# is always greater than or equal to min_delta. + +procedure ph_deltamin (a, min_delta, b, npix) + +real a[ARB] # input vector +real min_delta # minimum permitted absolute value +real b[ARB] # output vector +int npix # number of points + +int i + +begin + do i = 1, npix { + if (abs(a[i]) >= min_delta) + b[i] = a[i] + else if (a[i] < 0.0) + b[i] = -min_delta + else + b[i] = min_delta + } +end diff --git a/noao/digiphot/photcal/evaluate/phminv.f b/noao/digiphot/photcal/evaluate/phminv.f new file mode 100644 index 00000000..44d7430e --- /dev/null +++ b/noao/digiphot/photcal/evaluate/phminv.f @@ -0,0 +1,114 @@ +C SUBROUTINE PHMINV.F +C +C SOURCE +C BEVINGTON, PAGES 302-303. +C +C PURPOSE +C INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT +C +C USAGE +C CALL PHMINV (ARRAY, IK, JK, NORDER, DET) +C +C DESCRIPTION OF PARAMETERS +C ARRAY - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE +C IK - WORK ARRAY REQUIRED BY PHMINV +C JK - WORK ARRAY REQUIRED BY PHMINV +C NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT) +C DET - DETERMINANT OF INPUT MATRIX +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C NONE +C +C COMMENT +C DIMENSION STATEMENT VALID FOR NORDER UP TO 10 +C + SUBROUTINE PHMINV (ARRAY,IK,JK,NORDER,DET) + DOUBLE PRECISION ARRAY,AMAX,SAVE + DIMENSION ARRAY(NORDER,NORDER),IK(NORDER),JK(NORDER) +C +10 DET=1. +11 DO 100 K=1,NORDER +C +C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX +C + AMAX=0. +21 CONTINUE + IK(K) = K + JK(K) = K + DO 30 I=K,NORDER + DO 30 J=K,NORDER +23 IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30 +24 AMAX=ARRAY(I,J) + IK(K)=I + JK(K)=J +30 CONTINUE +C +C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) +C +41 I=IK(K) + IF (I-K) 21,51,43 +43 DO 50 J=1,NORDER + SAVE=ARRAY(K,J) + ARRAY(K,J)=ARRAY(I,J) +50 ARRAY(I,J)=-SAVE +51 J=JK(K) + IF (J-K) 21,600,53 +53 DO 60 I=1,NORDER + SAVE=ARRAY(I,K) + ARRAY(I,K)=ARRAY(I,J) +60 ARRAY(I,J)=-SAVE + +C +C CHECK FOR SINGULAR MATRIX +C +600 IF (AMAX) 61,610,61 +610 DET=0. + DO 700 I=1,NORDER + IF (I-K) 630,700,630 +630 ARRAY(I,K)=0. +700 CONTINUE + DO 800 J=1,NORDER + IF (J-K) 730,800,730 +730 ARRAY(K,J)=0. +800 CONTINUE + ARRAY(K,K)=0. + GOTO 100 +C +C ACCUMULATE ELEMENTS OF INVERSE MATRIX +C +61 DO 70 I=1,NORDER + IF (I-K) 63,70,63 +63 ARRAY(I,K)=-ARRAY(I,K)/AMAX +70 CONTINUE +71 DO 80 I=1,NORDER + DO 80 J=1,NORDER + IF (I-K) 74,80,74 +74 IF (J-K) 75,80,75 +75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J) +80 CONTINUE +81 DO 90 J=1,NORDER + IF (J-K) 83,90,83 +83 ARRAY(K,J)=ARRAY(K,J)/AMAX +90 CONTINUE + ARRAY(K,K)=1./AMAX +100 DET=DET*AMAX +C +C RESTORE ORDERING OF MATRIX +C +101 DO 130 L=1,NORDER + K=NORDER-L+1 + J=IK(K) + IF (J-K) 111,111,105 +105 DO 110 I=1,NORDER + SAVE=ARRAY(I,K) + ARRAY(I,K)=-ARRAY(I,J) +110 ARRAY(I,J)=SAVE +111 I=JK(K) + IF (I-K) 130,130,113 +113 DO 120 J=1,NORDER + SAVE=ARRAY(K,J) + ARRAY(K,J)=-ARRAY(I,J) +120 ARRAY(I,J)=SAVE +130 CONTINUE +140 RETURN + END diff --git a/noao/digiphot/photcal/evaluate/phprint.x b/noao/digiphot/photcal/evaluate/phprint.x new file mode 100644 index 00000000..94150792 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/phprint.x @@ -0,0 +1,110 @@ +include "../lib/parser.h" + +# PH_MKPLIST -- Construct the list of variables to be printed. + +int procedure ph_mkplist (plist, cmap, omap, nobsvars, psym, pcols, + max_len_plist) + +int plist # pointer to the list of variables +int cmap # catalog column to data column mapping +int omap # observations column to data column mapping +int nobsvars # number of observations variables +int psym[ARB] # the output array of variable symbols +int pcols[ARB] # offset in the data array of the symbol +int max_len_plist # the maximum length of the variables list + +int len_plist, sym, col +pointer sp, pname +int fntgfnb(), pr_getsym(), pr_gsymi(), pr_findmap1() + +begin + call smark (sp) + call salloc (pname, SZ_FNAME, TY_CHAR) + + len_plist = 0 + while (fntgfnb (plist, Memc[pname], SZ_FNAME) != EOF) { + if (len_plist >= max_len_plist) + break + sym = pr_getsym (Memc[pname]) + if (IS_INDEFI(sym)) + next + switch (pr_gsymi (sym, PSYMTYPE)) { + case PTY_CATVAR: + psym[len_plist+1] = sym + col = pr_gsymi (sym, PINPCOL) + pcols[len_plist+1] = pr_findmap1 (cmap, col) + nobsvars + case PTY_OBSVAR: + psym[len_plist+1] = sym + col = pr_gsymi (sym, PINPCOL) + pcols[len_plist+1] = pr_findmap1 (omap, col) + case PTY_SETEQ: + psym[len_plist+1] = sym + pcols[len_plist+1] = 0 + default: + next + } + len_plist = len_plist + 1 + } + + call sfree (sp) + return (len_plist) +end + + +# PH_OFORMAT -- Construct the output formatstr string. + +procedure ph_oformatstr (getid, ncols, formatstr, maxch) + +int getid # output the object id +int ncols # number of columns in the output text file +char formatstr[ARB] # the output formatstr string +int maxch # maximum number of characters in the formatstr string + +int i, fcol + +begin + if (getid == YES) { + call strcpy ("%-10s ", formatstr, maxch) + fcol = 2 + } else { + formatstr[1] = EOS + fcol = 1 + } + + do i = fcol, ncols { + if (i == ncols) + call strcat ("%-7.3f\n", formatstr, maxch) + else + call strcat ("%-7.3f ", formatstr, maxch) + } +end + + +# PH_OFIELDS -- Count the number of fields in the formatstr string. + +int procedure ph_ofields (formatstr) + +char formatstr[ARB] # the input formatstr string + +char percent +int ip, findex, nfields +int stridx() +data percent /'%'/ + +begin + nfields = 0 + + ip = 1 + while (formatstr[ip] != EOS) { + findex = stridx (percent, formatstr[ip]) + if (findex == 0) + break + ip = findex + ip + if (formatstr[ip] == percent) + ip = ip + 1 + else + nfields = nfields + 1 + } + + return (nfields) +end diff --git a/noao/digiphot/photcal/evaluate/t_evalfit.x b/noao/digiphot/photcal/evaluate/t_evalfit.x new file mode 100644 index 00000000..daf407c2 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/t_evalfit.x @@ -0,0 +1,495 @@ +include <error.h> +include <math/nlfit.h> +include "../lib/io.h" +include "../lib/parser.h" + +# Define the pointer Mem +define MEMP Memi + +# T_EVALFIT - Main photometry processing task. This task will convert +# intrumental photometric indices into standard indices, by using the +# same configuration table and coefficients found by FITCOEFFS. + +procedure t_evalfit() + +pointer observations # pointer to the observations files list +pointer catalogs # pointer to the catalogs files list +pointer config # pointer to configuration file name +pointer paramfile # pointer to fitted parameters file name +pointer calibfile # pointer to the output file name +int type # type of output to be processed +int etype # algorithm for computing errors +pointer print # pointer to the output variables list +pointer formatstr # pointer to the output format string +pointer catdir # pointer to the standard catalogs directory + +int i, j, getid, matchid, vcol, ecol, pindex, dummy, stat +int obslist, stdlist, plist, ofd, ifd, sym, symvar, ncols +int ntrneqs, nparams, nstd, nobsvars, nvars, len_plist +pointer sp, input, starname, dummyname, stdtable, cmap, omap +pointer vars, uservars, usererrs, fsym, rsym, esym, params, errors, psym, pcols +real ref, fit, errval, resid, chisqr, rms, pval + +bool clgetb() +int clgwrd(), open(), io_gcoeffs(), io_gobs(), pr_parse(), pr_geti() +int pr_gsym(), pr_gsymi(), pr_findmap1(), pr_gvari(), fntopnb() +int fntgfnb(), fntlenb(), ph_mkplist(), ph_header(), ph_ofields() +pointer pr_gsymp(), pr_xgetname() +real pr_eval(), ph_erval() +errchk io_gcoeffs() + +begin + # Allocate space for file names and character strings. + call smark (sp) + call salloc (observations, SZ_LINE, TY_CHAR) + call salloc (catalogs, SZ_LINE, TY_CHAR) + call salloc (config, SZ_FNAME, TY_CHAR) + call salloc (paramfile, SZ_FNAME, TY_CHAR) + call salloc (calibfile, SZ_FNAME, TY_CHAR) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (starname, SZ_LINE, TY_CHAR) + call salloc (dummyname, SZ_LINE, TY_CHAR) + call salloc (print, SZ_LINE, TY_CHAR) + call salloc (formatstr, SZ_LINE, TY_CHAR) + call salloc (catdir, SZ_FNAME, TY_CHAR) + + # Get the input and output file lists and names. + call clgstr ("observations", Memc[observations], SZ_LINE) + call clgstr ("config", Memc[config], SZ_FNAME) + call clgstr ("parameters", Memc[paramfile], SZ_FNAME) + call clgstr ("calib", Memc[calibfile], SZ_FNAME) + call clgstr ("catalogs", Memc[catalogs], SZ_LINE) + call clgstr ("catdir", Memc[catdir], SZ_FNAME) + + # Get output type flags. + etype = clgwrd ("errors", Memc[dummyname], SZ_LINE, ERR_OPTIONS) + type = clgwrd ("objects", Memc[dummyname], SZ_LINE, TYPE_STRING) + call clgstr ("print", Memc[print], SZ_LINE) + call clgstr ("format", Memc[formatstr], SZ_LINE) + + # Open the output file. + iferr { + if (clgetb ("append")) + ofd = open (Memc[calibfile], APPEND, TEXT_FILE) + else + ofd = open (Memc[calibfile], NEW_FILE, TEXT_FILE) + } then { + call erract (EA_WARN) + call sfree (sp) + return + } + + # Parse the configuration table. + if (pr_parse (Memc[config]) == ERR) { + call eprintf ("Error: Cannot parse the configuration file\n") + call close (ofd) + call sfree (sp) + return + } + + # Get number of variables and allocate memory for their values. + # Program stars do not have standard indices (last values in + # the array), and they will remain undefined, and this will give + # undefined equation values for equations that use them. + + nobsvars = pr_geti (NOBSVARS) + nvars = nobsvars + pr_geti (NCATVARS) + call salloc (vars, nvars, TY_REAL) + call salloc (uservars, nvars, TY_INT) + call aclri (Memi[uservars], nvars) + call salloc (usererrs, nvars, TY_INT) + call aclri (Memi[usererrs], nvars) + + # Map observational and catalog columns. This has to be done here in + # order to avoid mapping columns every time a new program + # star is read. + + call pr_catmap (cmap, dummy) + call pr_obsmap (omap, dummy) + + # Get all the parameter values for all the equations. This has + # to be done outside the main loop to minimize disk i/o. + # Otherwise it would be necessary to fetch them for each + # star observation. + + ntrneqs = pr_geti (NTRNEQS) + call salloc (fsym, ntrneqs, TY_POINTER) + call salloc (rsym, ntrneqs, TY_POINTER) + call salloc (esym, ntrneqs, TY_POINTER) + call salloc (params, ntrneqs, TY_POINTER) + call salloc (errors, ntrneqs, TY_POINTER) + do i = 1, ntrneqs { + + # Get equation symbol and number of parameters for it. + sym = pr_gsym (i, PTY_TRNEQ) + nparams = pr_gsymi (sym, PTEQNPAR) + + # Get the fit, reference and error equation pointers. + Memi[fsym+i-1] = pr_gsymp (sym, PTEQRPNFIT) + Memi[rsym+i-1] = pr_gsymp (sym, PTEQRPNREF) + Memi[esym+i-1] = pr_gsymp (sym, PTEQRPNERROR) + + # Determine which catalog and observational variables were + # actually used in the reference equations and which have + # defined error columns. + + do j = 1, pr_gsymi (sym, PTEQNRCAT) + pr_gsymi(sym, PTEQNROBS) { + symvar = pr_gvari (sym, j, PTEQREFVAR) + vcol = pr_gsymi (symvar, PINPCOL) + ecol = pr_gsymi (symvar, PINPERRCOL) + if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) { + vcol = pr_findmap1 (cmap, vcol) + nobsvars + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (cmap, ecol) + nobsvars + } else { + vcol = pr_findmap1 (omap, vcol) + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (omap, ecol) + } + Memi[uservars+vcol-1] = vcol + if (! IS_INDEFI(ecol)) + Memi[usererrs+vcol-1] = ecol + } + + # Determine which catalog and observational variables were + # actually used in the fit equations and which have + # defined error columns. + + do j = 1, pr_gsymi (sym, PTEQNFCAT) + pr_gsymi(sym, PTEQNFOBS) { + symvar = pr_gvari (sym, j, PTEQFITVAR) + vcol = pr_gsymi (symvar, PINPCOL) + ecol = pr_gsymi (symvar, PINPERRCOL) + if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) { + vcol = pr_findmap1 (cmap, vcol) + nobsvars + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (cmap, ecol) + nobsvars + } else { + vcol = pr_findmap1 (omap, vcol) + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (omap, ecol) + } + Memi[uservars+vcol-1] = vcol + if (! IS_INDEFI(ecol)) + Memi[usererrs+vcol-1] = ecol + } + + # Allocate space for parameter values and errors, + # for the current equation, and read them from the + # parameters file. + + call salloc (MEMP[params+i-1], nparams, TY_REAL) + call salloc (MEMP[errors+i-1], nparams, TY_REAL) + iferr { + if (io_gcoeffs (Memc[paramfile], sym, stat, chisqr, rms, + Memr[MEMP[params+i-1]], Memr[MEMP[errors+i-1]], + nparams) != nparams) { + call eprintf ("Warning: Error reading parameters for ") + call eprintf ("equation %s from %s\n") + call pargstr (Memc[pr_xgetname (sym)]) + call pargstr (Memc[paramfile]) + call amovkr (INDEFR, Memr[MEMP[params+i-1]], nparams) + call amovkr (INDEFR, Memr[MEMP[errors+i-1]], nparams) + } + } then { + call erract (EA_WARN) + call close (ofd) + call pr_unmap (cmap) + call pr_unmap (omap) + call pr_free() + call sfree (sp) + return + } + + # Issue a warning if some of the equations in the system did + # not converge but proceed with the evaluation. + # didn't converge. + + if (stat != DONE) { + call eprintf ( + "Warning: The solution for equation %s did not converge") + call pargstr (Memc[pr_xgetname (sym)]) + } + } + + # Decide whether to use id and/or catalog matching or not. + if (pr_geti (MINCOL) <= 1) { + getid = NO + matchid = NO + stdtable = NULL + } else if (Memc[catalogs] == EOS) { + getid = YES + matchid = NO + stdtable = NULL + } else { + getid = YES + matchid = YES + } + + # Get the list of optional variables to be printed. These variables + # may include any of the catalog or observations variables or the + # set equations. Variables which match any of these catagories + # will be ommited from the list to be printed. + + plist = fntopnb (Memc[print], NO) + len_plist = fntlenb (plist) + if (len_plist > 0) { + call salloc (psym, len_plist, TY_INT) + call salloc (pcols, len_plist, TY_INT) + len_plist = ph_mkplist (plist, cmap, omap, nobsvars, Memi[psym], + Memi[pcols], len_plist) + } + call fntclsb (plist) + + # Print the output header. + ncols = ph_header (ofd, Memc[observations], Memc[catalogs], + Memc[config], Memc[paramfile], type, getid, Memi[psym], + Memi[pcols], len_plist, etype, matchid) + + # Set the formatstr string. + if (Memc[formatstr] == EOS) { + call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE) + } else if (ph_ofields (Memc[formatstr]) != ncols) { + call eprintf ("Warning: The number of formatstr string fields ") + call eprintf ("does not match the number of output columns\n") + call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE) + } + + # If catalog matching is enabled, read in the catalog data. + if (matchid == YES) { + stdlist = fntopnb (Memc[catalogs], NO) + call io_gcatdat (Memc[catdir], stdlist, stdtable, nstd, dummy) + call fntclsb (stdlist) + } + + # Loop over the observation files. + obslist = fntopnb (Memc[observations], NO) + while (fntgfnb (obslist, Memc[input], SZ_FNAME) != EOF) { + + # Open the input file. + iferr (ifd = open (Memc[input], READ_ONLY, TEXT_FILE)) { + call erract (EA_WARN) + next + } + + # Get observations and names for program stars. The call + # to io_getline_init() is necessary since it's not called + # inside io_gobs(), and it should be called at least once + # per input file. + + call io_getline_init() + while (io_gobs (ifd, stdtable, omap, type, Memr[vars], nvars, + getid, Memc[starname], Memc[dummyname], SZ_LINE) != EOF) { + + # Print the star name. + call fprintf (ofd, Memc[formatstr]) + if (getid == YES) + call pargstr (Memc[starname]) + + # Print the extra variables. + do i = 1, len_plist { + pindex = Memi[pcols+i-1] + if (pindex > 0) + pval = Memr[vars+pindex-1] + else + pval = pr_eval (Memi[psym+i-1], Memr[vars], + Memr[MEMP[params]]) + call pargr (pval) + } + + # Loop over all the transformatstrion equations. + do i = 1, pr_geti (NTRNEQS) { + + # Get the equation symbol. + sym = pr_gsym (i, PTY_TRNEQ) + + # Compute the fit. + fit = pr_eval (Memi[fsym+i-1], Memr[vars], + Memr[MEMP[params+i-1]]) + call pargr (fit) + + # Compute the error. + if (etype != ERR_UNDEFINED) { + if (IS_INDEFR(fit)) + errval = INDEFR + switch (etype) { + case ERR_OBSERRORS: + if (IS_INDEFR(fit)) + errval = INDEFR + else + errval = ph_erval (Memi[fsym+i-1], fit, + Memr[MEMP[params+i-1]], Memr[vars], + Memi[uservars], Memi[usererrs], nobsvars) + case ERR_EQUATIONS: + if (Memi[esym+i-1] == NULL) + errval = INDEFR + else + errval = pr_eval (Memi[esym+i-1], Memr[vars], + Memr[MEMP[params+i-1]]) + default: + errval = INDEFR + } + call pargr (errval) + } + + # Compute the residual. + if (getid == NO || (matchid == YES && + type != TYPE_PROGRAM)) { + if (IS_INDEFR (fit)) + resid = INDEFR + else { + ref = pr_eval (Memi[rsym+i-1], Memr[vars], + Memr[MEMP[params+i-1]]) + if (IS_INDEFR(ref)) + resid = INDEFR + else + resid = ref - fit + } + call pargr (resid) + } + } + } + + # Close the input file. + call close (ifd) + } + + # Free memory. + call sfree (sp) + call pr_free() + call pr_unmap (cmap) + call pr_unmap (omap) + + # Close all files. + if (stdtable != NULL) + call stclose (stdtable) + call close (ofd) + call fntclsb (obslist) +end + + +# PH_HEADER - Print the output file header. + +int procedure ph_header (fd, observations, catalogs, config, paramfile, type, + getid, psym, pcols, len_plist, etype, matchid) + +int fd # output file descriptor +char observations[ARB] # observation file list +char catalogs[ARB] # catalog file list +char config[ARB] # configuration file name +char paramfile[ARB] # fitted parameters file +int type # type of object to output +int getid # output the object id +int psym[ARB] # list of additional variables to be output +int pcols[ARB] # columns numbers of additional variables +int len_plist # length of symbol list +int etype # type of error output +int matchid # is it possible to match ids + +int i, olist, slist, ncols +pointer sp, time, name +int fntopnb(), fntgfnb(), pr_geti(), pr_gsym() +long clktime() +pointer pr_gsymc(), pr_gsymp(), pr_xgetname() + +begin + # Allocate working space. + call smark (sp) + call salloc (time, SZ_LINE, TY_CHAR) + call salloc (name, SZ_FNAME, TY_CHAR) + + # Add the time stamp. + call cnvtime (clktime (0), Memc[time], SZ_LINE) + call fprintf (fd, "\n# %s\n") + call pargstr (Memc[time]) + + # Add the observation file names. + olist = fntopnb (observations, NO) + call fprintf (fd, "# List of observations files:\n") + while (fntgfnb (olist, Memc[name], SZ_FNAME) != EOF) { + call fprintf (fd, "#\t\t%s\n") + call pargstr (Memc[name]) + } + call fntclsb (olist) + + # Add the catalog file names. + if (matchid == YES) { + slist = fntopnb (catalogs, NO) + call fprintf (fd, "# List of catalog files:\n") + while (fntgfnb (slist, Memc[name], SZ_FNAME) != EOF) { + call fprintf (fd, "#\t\t%s\n") + call pargstr (Memc[name]) + } + call fntclsb (slist) + } + + # Add the configuration file name. + call fprintf (fd, "# Config:\t%s\n") + call pargstr (config) + + # Add the parameters file name. + call fprintf (fd, "# Parameters:\t%s\n") + call pargstr (paramfile) + + # Add the output options. + call fprintf (fd, "#\n") + if (matchid == YES) { + if (type == TYPE_ALL) + call fprintf (fd, + "# Computed indices for program and standard objects\n") + else if (type == TYPE_PROGRAM) + call fprintf (fd, + "# Computed indices for program objects only\n") + else if (type == TYPE_STANDARDS) + call fprintf (fd, + "# Computed indices for standard objects only\n") + } else + call fprintf (fd, + "# Computed indices for program and standard objects\n") + + # Print the optional id header. + call fprintf (fd, "#\n") + call fprintf (fd, "# Columns:\n") + if (getid == YES) { + call fprintf (fd, "#\t1\tobject id\n") + ncols = 1 + } else + ncols = 0 + + # Print headers for the variables in the print list and set any set + # equation offsets to pointers. + do i = 1, len_plist { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\t%s\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(psym[i])]) + if (pcols[i] <= 0) + psym[i] = pr_gsymp (psym[i], PSEQRPNEQ) + } + + + # Print the equation headers. + do i = 1, pr_geti (NTRNEQS) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\t%s\n") + call pargi (ncols) + call pargstr (Memc[pr_gsymc(pr_gsym (i,PTY_TRNEQ),PTEQREF)]) + if (etype != ERR_UNDEFINED) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\terror(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_gsymc(pr_gsym (i,PTY_TRNEQ),PTEQREF)]) + } + if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\tresid(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_gsymc (pr_gsym(i,PTY_TRNEQ),PTEQREF)]) + } + } + + call fprintf (fd,"\n\n") + + call sfree (sp) + + return (ncols) +end diff --git a/noao/digiphot/photcal/evaluate/t_invertfit.x b/noao/digiphot/photcal/evaluate/t_invertfit.x new file mode 100644 index 00000000..532d7f08 --- /dev/null +++ b/noao/digiphot/photcal/evaluate/t_invertfit.x @@ -0,0 +1,792 @@ +include <error.h> +include <math/nlfit.h> +include "../lib/io.h" +include "../lib/parser.h" +include "../lib/preval.h" + +# Define the pointer Mem +define MEMP Memi + +# T_INVERTFIT - INVERTFIT converts intrumental photometric indices into +# standard indices by using the configuration file, the coefficients +# determined by FITPARAMS and inverting the transformations. + +procedure t_invertfit () + +pointer observations # pointer to the observations file list +pointer catalogs # pointer to the catalogs file list +pointer config # pointer to configuration file name +pointer paramfile # pointer to fitted parameters file name +pointer calibfile # pointer to the output file name +int type # type of output to be processed +int etype # algorithm for computing the errors +pointer print # pointer to output variables list +pointer formatstr # pointer to the output format string +pointer catdir # pointer to the standard star directory + +int i, j, vcol, ecol, pindex, dummy, stat, getid, matchid +int obslist, stdlist, plist, ofd, ifd, sym, symvar, ncols, nstd, nset +int ntrneqs, nparams, nstdvars, nustdvars, nobsvars, nvars, nueq, maxnset +int len_plist, refcode +pointer sp, input, starname, dummyname, stdtable, omap, cmap +pointer vars, eqvartable, uservars, usererrs, userset, eqset, tvars +pointer dtvars, avtvars, ervars, svars, servars, params, errors, psym, pcols +pointer varindex, eqindex +real resid, chisqr, rms, pval + +bool clgetb() +int fntopnb(), fntgfnb(), clgwrd(), open(), io_gcoeffs(), io_gobs() +int pr_parse(), pr_geti(), pr_gsym(), pr_gsymi(), pr_gvari(), ph_objcheck() +int ph_invert(), pr_findmap1(), ph_iereqn(), ph_ierval(), ph_seteqn() +int ph_mkplist(), fntlenb(), ph_ofields(), ph_iheader(), ph_setvar() +pointer pr_xgetname(), pr_gsymp() +real pr_eval() +errchk io_gcoeffs() + +begin + # Allocate space for file names and character strings. + + call smark (sp) + call salloc (observations, SZ_LINE, TY_CHAR) + call salloc (catalogs, SZ_LINE, TY_CHAR) + call salloc (config, SZ_FNAME, TY_CHAR) + call salloc (paramfile, SZ_FNAME, TY_CHAR) + call salloc (calibfile, SZ_FNAME, TY_CHAR) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (starname, SZ_LINE, TY_CHAR) + call salloc (dummyname, SZ_LINE, TY_CHAR) + call salloc (print, SZ_LINE, TY_CHAR) + call salloc (formatstr, SZ_LINE, TY_CHAR) + call salloc (catdir, SZ_FNAME, TY_CHAR) + + # Get the observations list, the catalog list, the configuration + # file, the parameters file, and the output file names. + + call clgstr ("observations", Memc[observations], SZ_LINE) + call clgstr ("config", Memc[config], SZ_FNAME) + call clgstr ("parameters", Memc[paramfile], SZ_FNAME) + call clgstr ("calib", Memc[calibfile], SZ_FNAME) + call clgstr ("catalogs", Memc[catalogs], SZ_LINE) + call clgstr ("catdir", Memc[catdir], SZ_LINE) + + # Get the output type flags. + + etype = clgwrd ("errors", Memc[dummyname], SZ_LINE, ERR_OPTIONS) + type = clgwrd ("objects", Memc[dummyname], SZ_LINE, TYPE_STRING) + call clgstr ("print", Memc[print], SZ_LINE) + call clgstr ("format", Memc[formatstr], SZ_LINE) + + # Open the output file. + + iferr { + if (clgetb ("append")) + ofd = open (Memc[calibfile], APPEND, TEXT_FILE) + else + ofd = open (Memc[calibfile], NEW_FILE, TEXT_FILE) + } then { + call erract (EA_WARN) + call sfree (sp) + return + } + + # Parse the configuration file. + + if (pr_parse (Memc[config]) == ERR) { + call eprintf ("Error: Cannot parse the configuration file\n") + call close (ofd) + call sfree (sp) + return + } + + # Fetch the total number of observed and catalog variables and + # the total number of equations in the configuration file. + + nobsvars = pr_geti (NOBSVARS) + nstdvars = pr_geti (NCATVARS) + nvars = nobsvars + nstdvars + ntrneqs = pr_geti (NTRNEQS) + + # Map observations file and catalog file columns. + + call pr_catmap (cmap, dummy) + call pr_obsmap (omap, dummy) + + # Check which catalog and observations variables are actually + # used in the equations to be inverted and determine whether the + # system of equations can actually be inverted. This step is the + # first pass through the system of equations. A second pass is + # necessary to deal with the set equations. Use this first pass + # to fetch the fitted parameters for the equations. + + call salloc (eqvartable, nstdvars * ntrneqs, TY_INT) + call aclri (Memi[eqvartable], nstdvars * ntrneqs) + call salloc (uservars, nstdvars, TY_INT) + call aclri (Memi[uservars], nstdvars) + call salloc (usererrs, nobsvars, TY_INT) + call aclri (Memi[usererrs], nobsvars) + call salloc (params, ntrneqs, TY_POINTER) + call salloc (errors, ntrneqs, TY_POINTER) + + do i = 1, ntrneqs { + + # Get the equation symbol. + sym = pr_gsym (i, PTY_TRNEQ) + + # Get the reference equation symbol. + refcode = pr_gsymp (sym, PTEQRPNREF) + + # Quit if there are catalog variables in the function expression + # or if there are set equations containing references to the + # catalog variables in the reference equation. + + if ((pr_gsymi (sym, PTEQNRCAT) > 0) || (ph_seteqn (refcode) == + YES)) { + call eprintf ("Error: Cannot invert equations with catalog ") + call eprintf ("variables in the function expression\n") + call pr_unmap (cmap) + call pr_unmap (omap) + call close (ofd) + call sfree (sp) + return + } + + # Determine which catalog and observational variables were + # actually used in the reference equations and which have + # defined error columns. + + do j = 1, pr_gsymi (sym, PTEQNRCAT) + pr_gsymi(sym, PTEQNROBS) { + symvar = pr_gvari (sym, j, PTEQREFVAR) + vcol = pr_gsymi (symvar, PINPCOL) + ecol = pr_gsymi (symvar, PINPERRCOL) + if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) { + vcol = pr_findmap1 (cmap, vcol) + Memi[eqvartable+(i-1)*nstdvars+vcol-1] = vcol + nobsvars + Memi[uservars+vcol-1] = vcol + nobsvars + } else { + vcol = pr_findmap1 (omap, vcol) + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (omap, ecol) + if (! IS_INDEFI(ecol)) + Memi[usererrs+vcol-1] = ecol + } + } + + # Determine which catalog and observational variables were + # actually used in the fit equations and which have + # defined error columns. + + do j = 1, pr_gsymi (sym, PTEQNFCAT) + pr_gsymi(sym, PTEQNFOBS) { + symvar = pr_gvari (sym, j, PTEQFITVAR) + vcol = pr_gsymi (symvar, PINPCOL) + ecol = pr_gsymi (symvar, PINPERRCOL) + if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) { + vcol = pr_findmap1 (cmap, vcol) + Memi[eqvartable+(i-1)*nstdvars+vcol-1] = vcol + nobsvars + Memi[uservars+vcol-1] = vcol + nobsvars + } else { + vcol = pr_findmap1 (omap, vcol) + if (! IS_INDEFI(ecol)) + ecol = pr_findmap1 (omap, ecol) + if (! IS_INDEFI(ecol)) + Memi[usererrs+vcol-1] = ecol + } + } + + # Get the number of parameters for the equation. + # Allocate space for parameter values and errors, for the + # current equation, and read them from the coefficient file. + + nparams = pr_gsymi (sym, PTEQNPAR) + call salloc (MEMP[params+i-1], nparams, TY_REAL) + call salloc (MEMP[errors+i-1], nparams, TY_REAL) + iferr { + if (io_gcoeffs (Memc[paramfile], sym, stat, chisqr, rms, + Memr[MEMP[params+i-1]], Memr[MEMP[errors+i-1]], nparams) != + nparams) { + call eprintf ("Warning: Error reading parameters for ") + call eprintf ("equation %s from %s\n") + call pargstr (Memc[pr_xgetname(sym)]) + call pargstr (Memc[paramfile]) + call amovkr (INDEFR, Memr[MEMP[params+i-1]], nparams) + call amovkr (INDEFR, Memr[MEMP[errors+i-1]], nparams) + } + } then { + call erract (EA_WARN) + call close (ofd) + call pr_unmap (cmap) + call pr_unmap (omap) + call pr_free() + call sfree (sp) + return + } + + # Issue a warning if any of the equations in the system to be + # inverted did not converge but proceed with the fit. + + if (stat != DONE) { + call eprintf ( + "Warning: The solution for equation %s did not converge") + call pargstr (Memc[pr_xgetname(sym)]) + } + + } + + # Count the number of catalog variables used in the transformation + # equations. + + nustdvars = 0 + do i = 1, nstdvars { + if (Memi[uservars+i-1] != 0) + nustdvars = nustdvars + 1 + } + + # If the number of equations is less than the number of referenced + # catalog variables it is not possible to invert the transformation + # equations. + + if (nustdvars > ntrneqs) { + call eprintf ("Error: The number of equations is less than ") + call eprintf ("the number of unknowns\n") + call close (ofd) + call pr_unmap (cmap) + call pr_unmap (omap) + call pr_free() + call sfree (sp) + return + } + + # Loop over the transformation equations a second time searching for + # references to the set equations in the fit expressions. If a set + # equation reference is found, check to see whether it contains any + # reference to a catalog variable which is not referenced elsewhere + # in the transformation equations. Recompute the number of catalog + # variables used in the fit equations. + + maxnset = max (1, pr_geti (NSETEQS)) + call salloc (eqset, maxnset * ntrneqs, TY_INT) + call aclri (Memi[eqset], maxnset * ntrneqs) + nset = 0 + if (pr_geti (NSETEQS) > 0) { + + # Find the number of independent set equations. + call salloc (userset, pr_geti (NSETEQS), TY_INT) + call amovki (NULL, Memi[userset], pr_geti (NSETEQS)) + do i = 1, ntrneqs { + sym = pr_gsym (i, PTY_TRNEQ) + nset = nset + ph_setvar (i, sym, cmap, omap, Memi[eqvartable], + Memi[uservars], Memi[usererrs], nstdvars, nobsvars, + Memi[eqset], pr_geti (NSETEQS), Memi[userset], nset) + } + + # Is the system invertable? + if ((nustdvars + nset) > ntrneqs) { + call eprintf ("Error: The number of equations is less than ") + call eprintf ("the number of unknowns\n") + call close (ofd) + call pr_unmap (cmap) + call pr_unmap (omap) + call pr_free() + call sfree (sp) + return + } + + # Recompute the number of independent catalog variables. + nustdvars = 0 + do i = 1, nstdvars { + if (Memi[uservars+i-1] != 0) + nustdvars = nustdvars + 1 + } + + } + + + # Decide whether to use id/catalog matching or not. If id matching + # is enabled fetch the catalog data. + + if (pr_geti (MINCOL) <= 1) { + getid = NO + matchid = NO + stdtable = NULL + } else if (Memc[catalogs] == EOS) { + getid = YES + matchid = NO + stdtable = NULL + } else { + getid = YES + matchid = YES + stdtable = NULL + } + + # Get the list of optional variables to be printed. These variables + # may include any of the catalog or observations variables or the + # set equations. Variables which do not match any of these categories + # will be ommitted from the list to be printed. + + plist = fntopnb (Memc[print], NO) + len_plist = fntlenb (plist) + if (len_plist > 0) { + call salloc (psym, len_plist, TY_INT) + call salloc (pcols, len_plist, TY_INT) + len_plist = ph_mkplist (plist, cmap, omap, nobsvars, Memi[psym], + Memi[pcols], len_plist) + + } + call fntclsb (plist) + + # Print the output header. + ncols = ph_iheader (ofd, Memc[observations], Memc[catalogs], + Memc[config], Memc[paramfile], type, getid, Memi[psym], + Memi[pcols], len_plist, Memi[uservars], nstdvars, + Memi[userset], nset, etype, matchid) + + # Set the format string. + if (Memc[formatstr] == EOS) { + call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE) + } else if (ph_ofields (Memc[formatstr]) != ncols) { + call eprintf ("Warning: The number of format string fields ") + call eprintf ("does not match the number of output columns\n") + call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE) + } + + # Read in the catalog data. + if (matchid == YES) { + stdlist = fntopnb (Memc[catalogs], NO) + call io_gcatdat (Memc[catdir], stdlist, stdtable, nstd, dummy) + call fntclsb (stdlist) + } + + # Compute the initial values for the catalog variables and initial + # values for the delta-variable estimates. If a catalog is present + # the initial value of each catalog variable is the average value of + # that variable in the catalog. The standard deviation of the catalog + # variables is used as the initial value for the catalog variable + # increments. + + call salloc (avtvars, nstdvars, TY_REAL) + call salloc (dtvars, nstdvars, TY_REAL) + call ph_avstdvars (stdtable, Memi[uservars], Memr[avtvars], + Memr[dtvars], nstdvars) + + # Allocated space for the catalog variables and the fitted variables + # and their errors. + call salloc (vars, nvars, TY_REAL) + call salloc (tvars, nvars, TY_REAL) + if (nset > 0) + call salloc (svars, nset, TY_REAL) + call salloc (ervars, nstdvars, TY_REAL) + if (nset > 0) + call salloc (servars, nset, TY_REAL) + call salloc (varindex, nstdvars, TY_INT) + call salloc (eqindex, ntrneqs, TY_INT) + + # Initialize the fit inversion code. + call ph_ivinit (nstdvars, nustdvars, ntrneqs) + + # Loop over the observations files. + obslist = fntopnb (Memc[observations], NO) + while (fntgfnb (obslist, Memc[input], SZ_FNAME) != EOF) { + + # Open the input file. + iferr (ifd = open (Memc[input], READ_ONLY, TEXT_FILE)) { + call erract (EA_WARN) + next + } + + # Get observations and names for program stars. The call + # to io_getline_init() is necessary since it's not called + # inside io_gobs(), and it should be called at least once + # per input file. + + call io_getline_init() + while (io_gobs (ifd, stdtable, omap, type, Memr[vars], nvars, + getid, Memc[starname], Memc[dummyname], SZ_LINE) != EOF) { + + # Print star name. + call fprintf (ofd, Memc[formatstr]) + if (getid == YES) + call pargstr (Memc[starname]) + #call eprintf ("name=%s\n") + #call pargstr (Memc[starname]) + + # Print the extra variables. + do i = 1, len_plist { + pindex = Memi[pcols+i-1] + if (pindex > 0) + pval = Memr[vars+pindex-1] + else + pval = pr_eval (Memi[psym+i-1], Memr[vars], + Memr[MEMP[params]]) + call pargr (pval) + } + + # Initialize the fit. + call amovr (Memr[vars], Memr[tvars], nobsvars) + call amovr (Memr[avtvars], Memr[tvars+nobsvars], nstdvars) + + # Invert the transformations and compute the errors. + if (ph_objcheck (MEMP[params], Memr[tvars], Memi[eqvartable], + nstdvars, ntrneqs, Memi[eqset], maxnset, + Memi[varindex], nustdvars, Memi[eqindex], + nueq) == ERR) { + call amovkr (INDEFR, Memr[tvars+nobsvars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[svars], nset) + call amovkr (INDEFR, Memr[ervars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[servars], nset) + + } else if (ph_invert (MEMP[params], Memr[tvars], nobsvars, + Memr[dtvars], Memi[varindex], nstdvars, + nustdvars, Memi[eqindex], nueq) == ERR) { + + call amovkr (INDEFR, Memr[tvars+nobsvars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[svars], nset) + call amovkr (INDEFR, Memr[ervars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[servars], nset) + + } else { + + # Set any unfitted variables to INDEF. + do i = nobsvars + 1, nvars { + if (Memi[varindex+i-nobsvars-1] == 0) + Memr[tvars+i-1] = INDEFR + } + + # Evaluate the set equations. + if (nset > 0) { + do i = 1, nset + Memr[svars+i-1] = pr_eval (Memi[userset+i-1], + Memr[tvars], Memr[MEMP[params+i-1]]) + } + + # Evaluate the errors. + switch (etype) { + case ERR_UNDEFINED: + call amovkr (INDEFR, Memr[ervars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[servars], nset) + case ERR_EQUATIONS: + if (ph_iereqn (MEMP[params], Memr[tvars], nobsvars, + Memr[dtvars], Memi[varindex], Memr[ervars], + nstdvars, Memi[userset], Memr[svars], + Memr[servars], nset, nustdvars, Memi[eqindex], + nueq) <= 0) { + call amovkr (INDEFR, Memr[ervars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[servars], nset) + } + case ERR_OBSERRORS: + if (ph_ierval (MEMP[params], Memr[tvars], + Memi[usererrs], nobsvars, Memr[dtvars], + Memi[varindex], Memr[ervars], nstdvars, + Memi[userset], Memr[svars], Memr[servars], nset, + nustdvars, Memi[eqindex], nueq) <= 0) { + call amovkr (INDEFR, Memr[ervars], nstdvars) + if (nset > 0) + call amovkr (INDEFR, Memr[servars], nset) + } + default: + call amovkr (INDEFR, Memr[ervars], nstdvars) + } + } + + # Write out the standard indices, errors and residuals. + do i = nobsvars + 1, nvars { + + # Skip catalog variables not used in the equation. + if (Memi[uservars+i-nobsvars-1] <= 0) + next + + call pargr (Memr[tvars+i-1]) + if (etype != ERR_UNDEFINED) { + if (IS_INDEFR(Memr[tvars+i-1])) + call pargr (INDEFR) + else + call pargr (Memr[ervars+i-nobsvars-1]) + } + + # Compute the residual of the fit. + if (getid == NO || (matchid == YES && + type != TYPE_PROGRAM)) { + if (IS_INDEFR (Memr[vars+i-1]) || + IS_INDEFR (Memr[tvars+i-1])) + resid = INDEFR + else + resid = Memr[vars+i-1] - Memr[tvars+i-1] + call pargr (resid) + } + } + + # Write out the set equations. + do i = 1, nset { + + # Write out the set equation. + call pargr (Memr[svars+i-1]) + + # Evaluate the error. This is tricky at the moment. + if (etype != ERR_UNDEFINED) { + if (IS_INDEFR(Memr[svars+i-1])) + call pargr (INDEFR) + else + call pargr (Memr[servars+i-1]) + } + + # Compute the residual of the fit. + if (getid == NO || (matchid == YES && + type != TYPE_PROGRAM)) { + resid = pr_eval (Memi[userset+i-1], Memr[vars], + Memr[MEMP[params+i-1]]) + if (IS_INDEFR (Memr[svars+i-1]) || IS_INDEFR (resid)) + resid = INDEFR + else + resid = resid - Memr[svars+i-1] + call pargr (resid) + } + } + } + + # Close the input file. + call close (ifd) + } + + # Free memory. + call sfree (sp) + call ph_ivfree() + call pr_free() + call pr_unmap (cmap) + call pr_unmap (omap) + + # Close all the files. + if (stdtable != NULL) + call stclose (stdtable) + call close (ofd) + call fntclsb (obslist) +end + + +# PH_IHEADER - Print the output file header. + +int procedure ph_iheader (fd, observations, catalogs, config, paramfile, type, + getid, psym, pcols, len_plist, catvars, ncatvars, userset, nset, + etype, matchid) + +int fd # output file descriptor +char observations[ARB] # the list of observations files +char catalogs[ARB] # the list of catalog files +char config[ARB] # the configuration file name +char paramfile[ARB] # the fitted parameters file +int type # the type of object to output +int getid # output the object id ? +int psym[ARB] # list of additional symbols to be printed +int pcols[ARB] # column numbers of additional output variables +int len_plist # length of symbol list +int catvars[ARB] # the list of active catalog variables +int ncatvars # number of catalog variables +int userset[ARB] # list of set equations included in fit +int nset # number of set equations +int etype # type of error output +int matchid # is it possible to match ids ? + +int i, olist, slist, ncols +pointer sp, time, name, sym +int fntopnb(), fntgfnb(), pr_gsym() +long clktime() +pointer pr_xgetname(), pr_gsymp() + +begin + # Allocate some working space. + call smark (sp) + call salloc (time, SZ_LINE, TY_CHAR) + call salloc (name, SZ_FNAME, TY_CHAR) + + # Add the time stamp. + call cnvtime (clktime (0), Memc[time], SZ_LINE) + call fprintf (fd, "\n# %s\n") + call pargstr (Memc[time]) + + # Add the observation file names. + olist = fntopnb (observations, NO) + call fprintf (fd, "# List of observations files:\n") + while (fntgfnb (olist, Memc[name], SZ_FNAME) != EOF) { + call fprintf (fd, "#\t\t%s\n") + call pargstr (Memc[name]) + } + call fntclsb (olist) + + # Add the catalog file names. + if (matchid == YES) { + slist = fntopnb (catalogs, NO) + call fprintf (fd, "# Number of catalog files:\n") + while (fntgfnb (slist, Memc[name], SZ_FNAME) != EOF) { + call fprintf (fd, "#\t\t%s\n") + call pargstr (Memc[name]) + } + call fntclsb (slist) + } + + # Add the configuration file name. + call fprintf (fd, "# Config:\t%s\n") + call pargstr (config) + + # Add the parameters file name. + call fprintf (fd, "# Parameters:\t%s\n") + call pargstr (paramfile) + + # Write the output options. + call fprintf (fd, "#\n") + if (matchid == YES) { + if (type == TYPE_ALL) + call fprintf (fd, + "# Computed indices for program and standard objects\n") + else if (type == TYPE_PROGRAM) + call fprintf (fd, + "# Computed indices for program objects only\n") + else if (type == TYPE_STANDARDS) + call fprintf (fd, + "# Computed indices for standard objects only\n") + } else + call fprintf (fd, + "# Computed indices for program and standard objects\n") + + # Print the optional id header + call fprintf (fd, "#\n") + call fprintf (fd, "# Columns: \n") + if (getid == YES) { + call fprintf (fd, "#\t1\tobject id\n") + ncols = 1 + } else + ncols = 0 + + # Print headers for the variables in the print list and set any set + # equation offsets to pointers. + do i = 1, len_plist { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\t%s\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(psym[i])]) + if (pcols[i] <= 0) + psym[i] = pr_gsymp (psym[i], PSEQRPNEQ) + } + + # Print the catalog variables headers. + do i = 1, ncatvars { + if (catvars[i] <= 0) + next + sym = pr_gsym (i, PTY_CATVAR) + ncols = ncols + 1 + call fprintf (fd, "#\t%d\t%s\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(sym)]) + if (etype != ERR_UNDEFINED) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\terror(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(sym)]) + } + if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\tresid(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(sym)]) + } + } + + # Print the set equation variables headers. + do i = 1, nset { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\t%s\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(userset[i])]) + if (etype != ERR_UNDEFINED) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\terror(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(userset[i])]) + } + if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) { + ncols = ncols + 1 + call fprintf (fd, "#\t%d\tresid(%s)\n") + call pargi (ncols) + call pargstr (Memc[pr_xgetname(userset[i])]) + } + userset[i] = pr_gsymp (userset[i], PSEQRPNEQ) + } + + call fprintf (fd, "\n\n") + + call sfree (sp) + + return (ncols) +end + + +# PH_AVSTDVARS -- Compute the mean and sigma of the variables in the +# standard table to use as initial guesses for the math routines. +# If the standard catalog is undefined or empty the mean values of the +# variables are set to 0.0 and the sigmas are set to 1.0. If there is only +# one catalog value the sigma is also set to 1.0. + +procedure ph_avstdvars (stable, index, avg, sigma, npts) + +pointer stable # pointer to the symbol table +int index[ARB] # index of active catalog variables +real avg[ARB] # the output array of averages +real sigma[ARB] # the output array of standard deviations +int npts # number of points + +int i, ndata +pointer sym, sp, n +real data +pointer sthead(), stnext() + +begin + # Initialize and return if the standard table is undefined. + call amovkr (0.0, avg, npts) + if (stable == NULL) { + call amovkr (0.1, sigma, npts) + return + } + + # Allocate some temporary space and intialize the accumulators. + call smark (sp) + call salloc (n, npts, TY_INT) + call aclri (Memi[n], npts) + call aclrr (sigma, npts) + + # Read the symbol table and accumlate the sums. + sym = sthead (stable) + while (sym != NULL) { + do i = 1, npts { + if (index[i] == 0) + next + data = Memr[P2R(sym+i-1)] + if (IS_INDEFR(data)) + next + avg[i] = avg[i] + data + sigma[i] = sigma[i] + data ** 2 + Memi[n+i-1] = Memi[n+i-1] + 1 + } + sym = stnext (stable, sym) + } + + # Compute the averages. + do i = 1, npts { + ndata = Memi[n+i-1] + if (ndata == 0) + sigma[i] = 0.1 + else if (ndata == 1) + sigma[i] = 0.1 + else { + avg[i] = avg[i] / ndata + sigma[i] = (sigma[i] - avg[i] * avg[i] * ndata) / (ndata - 1) + if (sigma[i] <= 0.0) + sigma[i] = 0.1 + else + sigma[i] = min (0.1, sqrt (sigma[i])) + } + } + + call sfree (sp) +end |