aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/photcal/fitparams/fttrneq.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/photcal/fitparams/fttrneq.x')
-rw-r--r--noao/digiphot/photcal/fitparams/fttrneq.x541
1 files changed, 541 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/fitparams/fttrneq.x b/noao/digiphot/photcal/fitparams/fttrneq.x
new file mode 100644
index 00000000..6cfb092a
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/fttrneq.x
@@ -0,0 +1,541 @@
+include <error.h>
+include <pkg/xtanswer.h>
+include <pkg/gtools.h>
+include <math/nlfit.h>
+include <pkg/inlfit.h>
+include "../lib/parser.h"
+include "../lib/fitparams.h"
+
+# Interactive commands
+
+define CMD_OPTIONS "|first|last|prev|next|again|go|quit|"
+define CMD_FIRST 1
+define CMD_LAST 2
+define CMD_PREV 3
+define CMD_NEXT 4
+define CMD_AGAIN 5
+define CMD_GO 6
+define CMD_QUIT 7
+define CMD_MIN CMD_FIRST
+define CMD_MAX CMD_QUIT
+
+# Prompt to save fit results into output file and equation structure
+
+define SAVE_PROMPT "Do you want to save fit results ?"
+define QUIT_PROMPT "Program will quit. Do you really want to do so ?"
+
+# Program labels
+
+define retry 9999
+
+
+# FT_TRNEQS - Fit all the transformation equations.
+
+procedure ft_trneqs (output, logfile, graphics, otable, ntable, wtflag,
+ addscatter, tol, itmax, interactive, high, low, niterate, grow,
+ log_fit, log_results)
+
+char output[ARB] # output file name
+char logfile[ARB] # the log file name
+char graphics[ARB] # graphics output device
+pointer otable # standard observation table
+pointer ntable # standard name table
+int wtflag # the type of weighting
+int addscatter # add a scatter term to the weights
+real tol # fit tolerance
+int itmax # max number of iterations
+bool interactive # interactive fit ?
+real high, low # rejection factors
+int niterate # number of rejection iterations
+real grow # rejection growing radius
+bool log_fit # log the fit statistics
+bool log_results # log the fit results
+
+int cmd, neq, maxeq, nobs, sym
+pointer rtable, wtable, totable, gp, gt
+
+#bool clgetb()
+int mct_nrows(), pr_gsym(), pr_geti(), ft_trnanswer()
+pointer gopen(), gt_init(), pr_gsymp()
+
+begin
+ # Debug ?
+ #if (clgetb ("debug.fitcode")) {
+ #call eprintf ("ft_trneqs (output=%s) (graphics=%s) (ot=%d) "
+ #call pargstr (output)
+ #call pargstr (graphics)
+ #call pargi (otable)
+ #call eprintf ("(tol=%g) (maxiter=%d) (int=%b)\n")
+ #call pargr (tol)
+ #call pargi (itmax)
+ #call pargb (interactive)
+ #call eprintf ("ft_trneqs (gp=%d) (gt=%d) (high=%g) (lo=%g) ")
+ #call pargi (gp)
+ #call pargi (gt)
+ #call pargr (high)
+ #call pargr (low)
+ #call eprintf ("(niter=%d) (grow=%g)\n")
+ #call pargi (niterate)
+ #call pargr (grow)
+ #}
+
+ # Get number of observations.
+ nobs = mct_nrows (otable)
+
+ # Allocate space for reference and weight tables. These tables
+ # store the values of the reference equation for all observations,
+ # and the weight for each of them, respectively.
+
+ totable = NULL
+ call mct_alloc (rtable, nobs, 1, TY_REAL)
+ call mct_alloc (wtable, nobs, 1, TY_REAL)
+
+ # Open graphics.
+
+ if (interactive) {
+ gp = gopen (graphics, NEW_FILE, STDGRAPH)
+ gt = gt_init()
+ } else {
+ gp = NULL
+ gt = NULL
+ }
+
+ # Loop over transformation equations, until the user
+ # sends a quit command.
+
+ neq = 1
+ cmd = CMD_NEXT
+ maxeq = pr_geti (NTRNEQS)
+ repeat {
+
+ # Get the equation symbol.
+
+ sym = pr_gsym (neq, PTY_TRNEQ)
+
+ # Evaluate the reference equation for catalog
+ # observations and store them into a table.
+
+ call ft_rfeval (pr_gsymp (sym, PTEQRPNREF), otable, rtable)
+
+ # Evaluate weights for the current equation. If an error
+ # condition is raised (due to a negative weight value) the
+ # program continues with the next equation in the list.
+
+ switch (wtflag) {
+ case FWT_UNIFORM:
+ call ft_wtuniform (otable, wtable)
+
+ case FWT_PHOTOMETRIC:
+ iferr (call ft_wtphoterrs (sym, otable, wtable)) {
+ call erract (EA_WARN)
+ next
+ }
+ case FWT_EQUATIONS:
+ iferr (call ft_wteqns (sym, otable, wtable)) {
+ call erract (EA_WARN)
+ next
+ }
+ default:
+ call ft_wtuniform (otable, wtable)
+ }
+
+ # Save original observations, and check for undefined values.
+
+ call mct_copy (otable, totable)
+ call ft_indef (sym, totable, rtable, wtable)
+
+ # Fit transformation equation.
+
+ call ft_trneq (output, logfile, sym, totable, wtable, rtable,
+ ntable, wtflag, addscatter, tol, itmax, interactive, gp, gt,
+ high, low, niterate, grow, log_fit, log_results)
+
+ # Reset the reference table.
+
+ call mct_reset (rtable)
+
+retry
+ # Prompt the user for command.
+ if (interactive)
+ call ft_trncmd (cmd)
+ else
+ cmd = CMD_NEXT
+
+ # Branch on command.
+ switch (cmd) {
+ case CMD_FIRST:
+ neq = 1
+ case CMD_LAST:
+ neq = maxeq
+ case CMD_PREV:
+ if (neq > 1)
+ neq = neq - 1
+ else
+ goto retry
+ case CMD_NEXT:
+ neq = neq + 1
+ case CMD_AGAIN:
+ ;
+ case CMD_GO:
+ interactive = false
+ neq = neq + 1
+ case CMD_QUIT:
+ if (ft_trnanswer (QUIT_PROMPT, YES) == YES)
+ break
+ default:
+ call error (0, "fttrneqs: Unknown command")
+ }
+
+ # Prompt to quit.
+ if (neq > maxeq) {
+ if (interactive && cmd != CMD_GO) {
+ if (ft_trnanswer (QUIT_PROMPT, YES) == NO) {
+ neq = maxeq
+ goto retry
+ } else
+ break
+ } else
+ break
+ }
+ }
+
+ # Close graphics.
+ if (gp != NULL)
+ call gclose (gp)
+ if (gt != NULL)
+ call gt_free (gt)
+
+ # Free tables.
+ call mct_free (rtable)
+ call mct_free (wtable)
+ call mct_free (totable)
+end
+
+
+# FT_TRNEQ - Fit single transformation equation.
+
+procedure ft_trneq (output, logfile, sym, otable, wtable, rtable, ntable,
+ wtflag, addscatter, tol, itmax, interactive, gp, gt high, low,
+ niterate, grow, log_fit, log_results)
+
+char output[ARB] # output file name
+char logfile[ARB] # log file name
+int sym # transformation equation symbol
+pointer otable # standard observation table
+pointer wtable # weight table
+pointer rtable # reference equation table
+pointer ntable # names table
+int wtflag # type of weighting
+int addscatter # add a scatter term to the weight equation
+real tol # fit tolerance
+int itmax # max number of iterations
+bool interactive # interactive fit ?
+pointer gp, gt # GIO and GTOOLS descriptors
+real high, low # rejection factors
+int niterate # number of rejection iterations
+real grow # rejection growing radius
+bool log_fit # log the fit statistics
+bool log_results # log the fit results
+
+int nobs, nvars, nparams, nfparams, len_name
+int i, stat, answer, psym, nlwtflag
+pointer sp, params, dparams, name, aux, plist, nl, in
+real chisqr, variance, scatter, rms
+
+#bool clgetb()
+bool streq()
+int locpr(), mct_nrows(), mct_maxcol(), pr_geti(), pr_gsym(), pr_gsymi()
+int pr_gpari(), ft_trnanswer()
+pointer pr_xgetname(), pr_gsymp(), pr_gsymc(), mct_getbuf(), in_getp()
+real pr_gsymr()
+extern ft_func(), ft_dfunc(), ft_plot()
+
+begin
+ # Debug ?
+ #if (clgetb ("debug.fitcode")) {
+ #call eprintf ("ft_trneq (output=%s) (sym=%d) (ot=%d) (wt=%d) "
+ #call pargstr (output)
+ #call pargi (sym)
+ #call pargi (otable)
+ #call pargi (wtable)
+ #call eprintf ("(rt=%d) (tol=%g) (maxiter=%d)\n")
+ #call pargi (rtable)
+ #call pargr (tol)
+ #call pargi (itmax)
+ #call eprintf ("ft_trneq (int=%b) (gp=%d) (gt=%d) (high=%g) "
+ #call pargb (interactive)
+ #call pargi (gp)
+ #call pargi (gt)
+ #call pargr (high)
+ #call eprintf ("(low=%g) (niter=%d) (grow=%g)\n")
+ #call pargr (low)
+ #call pargi (niterate)
+ #call pargr (grow)
+ #}
+
+ # Get number of observations, and variables for the observational data.
+
+ nobs = mct_nrows (otable)
+ nvars = mct_maxcol (otable)
+ len_name = mct_maxcol (ntable)
+
+ # Get number of parameters for the current equation.
+
+ nparams = pr_gsymi (sym, PTEQNPAR)
+
+ # Allocate stack space.
+ call smark (sp)
+ call salloc (aux, SZ_LINE, TY_CHAR)
+ call salloc (name, SZ_LINE, TY_CHAR)
+ call salloc (params, nparams, TY_REAL)
+ call salloc (dparams, nparams, TY_REAL)
+ call salloc (plist, nparams, TY_INT)
+
+ # Initialize parameters, errors, and parameter list.
+
+ call aclrr (Memr[params], nparams)
+ call aclrr (Memr[dparams], nparams)
+ do i = 1, nparams {
+ psym = pr_gpari (sym, i, PTEQPAR)
+ Memr[dparams+i-1] = pr_gsymr (psym, PFITDELTA)
+ }
+ call amovi (Memi[pr_gsymp (sym, PTEQSPLIST)], Memi[plist], nparams)
+
+ # Initialize the fit evaluation. This is necessary in order to
+ # set the equations called by the INLFIT procedures.
+
+ call ft_evinit (sym, nparams)
+
+ # Get number of fitting parameters for the current equation.
+ nfparams = pr_gsymi (sym, PTEQNFPAR)
+
+ # Initialize INLFIT.
+ call in_initr (in, locpr (ft_func), locpr (ft_dfunc),
+ Memr[pr_gsymp (sym, PTEQSPARVAL)], Memr[dparams],
+ nparams, Memi[plist], nfparams)
+
+ # Set INLFIT fitting parameters.
+ call in_putr (in, INLTOLERANCE, tol)
+ call in_puti (in, INLMAXITER, itmax)
+ call in_putr (in, INLHIGH, high)
+ call in_putr (in, INLLOW, low)
+ call in_puti (in, INLNREJECT, niterate)
+ call in_putr (in, INLGROW, grow)
+
+ # Put in the reference and fit equation names.
+ call sprintf (Memc[aux], SZ_LINE, "|%s|%s|")
+ call pargstr (Memc[pr_gsymc (sym, PTEQREF)])
+ call pargstr (Memc[pr_gsymc (sym, PTEQFIT)])
+ call in_pstr (in, INLFLABELS, Memc[aux])
+
+ # Put in the parameter names.
+ call strcpy ("|", Memc[aux], SZ_LINE)
+ do i = 1, nparams {
+ call strcat (Memc[pr_xgetname (pr_gpari (sym, i, PTEQPAR))],
+ Memc[aux], SZ_LINE)
+ call strcat ("|", Memc[aux], SZ_LINE)
+ }
+ call in_pstr (in, INLPLABELS, Memc[aux])
+
+ # Put in the variable names.
+ call strcpy ("|", Memc[aux], SZ_LINE)
+ do i = 1, pr_geti (NOBSVARS) {
+ call pr_vtran (Memc[pr_xgetname (pr_gsym (i, PTY_OBSVAR))],
+ Memc[name], SZ_LINE)
+ call strcat (Memc[name], Memc[aux], SZ_LINE)
+ call strcat ("|", Memc[aux], SZ_LINE)
+ }
+ do i = 1, pr_geti (NCATVARS) {
+ call pr_vtran (Memc[pr_xgetname (pr_gsym (i, PTY_CATVAR))],
+ Memc[name], SZ_LINE)
+ call strcat (Memc[name], Memc[aux], SZ_LINE)
+ call strcat ("|", Memc[aux], SZ_LINE)
+ }
+ call in_pstr (in, INLVLABELS, Memc[aux])
+
+ # Put plot equations and redefine graph keys, but only if
+ # both plot equations are defined. Otherwise leave the defaults.
+
+ if (pr_gsymp (sym, PTEQRPNXPLOT) != NULL &&
+ pr_gsymp (sym, PTEQRPNYPLOT) != NULL) {
+
+ # Put in the plot equation names.
+ call sprintf (Memc[aux], SZ_LINE, "|%s|%s|")
+ call pargstr (Memc[pr_gsymc (sym, PTEQXPLOT)])
+ call pargstr (Memc[pr_gsymc (sym, PTEQYPLOT)])
+ call in_pstr (in, INLUSERLABELS, Memc[aux])
+
+ # Put in the plot equation.
+ call in_puti (in, INLUAXES, locpr (ft_plot))
+
+ # Redefine the last graph key.
+ call in_puti (in, INLGKEY, INLNGKEYS)
+ call in_pkey (in, INLNGKEYS, 1, KEY_UAXIS, 1)
+ call in_pkey (in, INLNGKEYS, 2, KEY_UAXIS, 2)
+ }
+
+ # Call the appropiate version of fitting routine according to the
+ # value of the interactive flag.
+
+ if (wtflag == FWT_UNIFORM)
+ nlwtflag = WTS_USER
+ else if (addscatter == YES)
+ nlwtflag = WTS_SCATTER
+ else
+ nlwtflag = WTS_USER
+
+ if (interactive) {
+ call gt_setr (gt, GTXMIN, INDEFR)
+ call gt_setr (gt, GTXMAX, INDEFR)
+ call gt_setr (gt, GTYMIN, INDEFR)
+ call gt_setr (gt, GTYMAX, INDEFR)
+ call ing_fitr (in, gp, "cursor", gt, nl, Memr[mct_getbuf (otable)],
+ Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)],
+ Memc[mct_getbuf(ntable)], nobs, nvars, len_name, nlwtflag,
+ stat)
+ } else {
+ call in_fitr (in, nl, Memr[mct_getbuf (otable)],
+ Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)],
+ nobs, nvars, nlwtflag, stat)
+ }
+
+ # Decide whether to prompt the user for saving answer.
+ if (interactive)
+ answer = ft_trnanswer (SAVE_PROMPT, YES)
+ else
+ answer = YES
+
+ # Write the parameter values into output file.
+ if (answer == YES || answer == ALWAYSYES) {
+
+ # Get parameter values, and parameter list.
+ call nlpgetr (nl, Memr[params], nparams)
+ call amovi (Memi[in_getp (in, INLPLIST)], Memi[plist], nparams)
+
+ # Get the parameter errors.
+ call in_errorsr (in, nl, Memr[mct_getbuf (otable)],
+ Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)], nobs,
+ nvars, variance, chisqr, scatter, rms, Memr[dparams])
+
+ # Write parameters and errors into output file.
+ iferr (call io_pcoeffs (output, sym, stat, wtflag, variance, chisqr,
+ scatter, rms, Memr[params], Memr[dparams], Memi[plist],
+ nparams))
+ call erract (EA_WARN)
+
+ # Log the fit and results.
+ if (log_fit || log_results) {
+ if (interactive && streq (logfile, "STDOUT"))
+ call gdeactivate (gp, 0)
+ if (log_fit) {
+ iferr (call io_title (logfile, "#EQUATION:", sym))
+ ;
+ iferr (call ing_showr (in, logfile))
+ ;
+ iferr (call ing_errorsr (in, logfile, nl,
+ Memr[mct_getbuf(otable)], Memr[mct_getbuf(rtable)],
+ Memr[mct_getbuf(wtable)], nobs, nvars))
+ ;
+ }
+ if (log_results) {
+ iferr (call io_title (logfile, "#RESULTS:", sym))
+ ;
+ iferr (call ing_resultsr (in, logfile, nl,
+ Memr[mct_getbuf(otable)], Memr[mct_getbuf(rtable)],
+ Memr[mct_getbuf(wtable)], Memc[mct_getbuf(ntable)],
+ nobs, nvars, len_name))
+ ;
+ }
+ if (interactive && streq (logfile, "STDOUT"))
+ call greactivate (gp, 0)
+ }
+
+
+ # Update fitted parameter values into the equation symbol
+ # substructure if the fit was succesfull.
+ if (stat == DONE)
+ call amovr (Memr[params],
+ Memr[P2R(pr_gsymp (sym, PTEQSPARVAL))], nparams)
+ }
+
+ # Debug ?
+ #if (clgetb ("debug.nlfit"))
+ #call dg_inldump (in, nl)
+
+
+ # Free inlfit and nlfit descriptors.
+ call in_freer (in)
+ call nlfreer (nl)
+
+ # Free fit evaluation.
+ call ft_evfree ()
+
+ # Free stack space.
+ call sfree (sp)
+end
+
+
+# FT_TRNCMD -- Prompt the user for a command.
+
+procedure ft_trncmd (cmd)
+
+int cmd # command code
+
+char command[SZ_LINE]
+int scan(), strlen()
+int strdic(), io_strwrd()
+
+begin
+ # Get current command string.
+ if (io_strwrd (cmd, command, SZ_LINE, CMD_OPTIONS) == 0)
+ call error (0, "ft_trneqs: Unknown command")
+
+ # Keep prompting the user until there is a valid command.
+ repeat {
+
+ # Print current setting.
+ call printf ("Command (%s) (%s) : ")
+ call pargstr (CMD_OPTIONS)
+ call pargstr (command)
+ call flush (STDOUT)
+
+ # Get new command.
+ if (scan () == EOF)
+ command[1] = EOS
+ else
+ call gargstr (command, SZ_LINE)
+ if (strlen (command) == 0) {
+ if (io_strwrd (cmd, command, SZ_LINE, CMD_OPTIONS) <= 0)
+ ;
+ }
+
+ # Search for command in the dictionary.
+ cmd = strdic (command, command, SZ_LINE, CMD_OPTIONS)
+
+ # Test command.
+ if (cmd >= CMD_MIN && cmd <= CMD_MAX)
+ break
+ else
+ call printf ("\007")
+ }
+end
+
+
+# FT_TRNANSWER -- Get a YES/NO answer from the user, and return it
+# as the procedure value.
+
+int procedure ft_trnanswer (prompt, default)
+
+char prompt # prompt to the user
+int default # default answer
+
+int answer
+
+begin
+ answer = default
+ call xt_answer (prompt, answer)
+ if (answer == YES || answer == ALWAYSYES)
+ return (YES)
+ else
+ return (NO)
+end
generated by cgit (git 2.34.1) at 2025-09-20 13:48:45 -0400