#@(#) t_autospec.x 2.18
include	<error.h>
include "../aspect.h"
include "../as_comp.h"
include "../as_err.h"
include "../as_fit.h"
include	"../as_spec.h"
include	"../as_units.h"

define 	DEBUG		false
define	MAX_DOMAINS	64

#--------------------------------------------------------------------10 Mar 96--
.help t_autospec.x Jan96 source/tasks
.ih
NAME
.nf
    t_autospec - Run the ASpect fitting in batch mode
get_cl_from_db - Build input component list from disk file & append to input CL
    fit_report - Report fit progress to STDOUT
.fi
.ih 
DESCRIPTION
Task to fit one or more model components to one or more input spectra. 
Fetches task parameters, initializes data structures, etc.
.ih
ROUTINE DETAILS
.ls procedure t_autospec ()
.ls ARGUMENTS
.le
.le

.ls procedure get_cl_from_db (cl, db_table)
.ls ARGUMENTS
.ls cl (input/output: pointer)
component list
.le
.ls db_table[ARB] (input: char)
name of input DB table of components
.le
.le
.le

.ls procedure fit_report (iter, chisq, force_flush)
.ls ARGUMENTS
.ls iter (input: int)
the current iteration count
.le
.ls chisq (input: real)
the current chisq value
.le
.ls force_flush (input: bool)
force flushing the output?
.le
.le
.le
.endhelp
#-------------------------------------------------------------------------------
#  T_AUTOSPEC -	Test driver for spectrum & component input from disk tables, 
#		and fitting routines.  Fetches task parameters, initializes 
#		data structures, etc.  
#
#  Revision history:
#	 2-Aug-94 by RAShaw	Initial implementation
#	 3-Aug-94 by RAShaw	Allow for multiple fit intervals
#	11-Aug-94 by RAShaw	Revise calling sequence for DB_READ to take 
#				a pointer to component array, and to return 
#				the no. component objects created; also, add 
#				error checking.  
#	12-Aug-94 by RAShaw	Move "fit" related structure from "aspect.h" to 
#				"as_fit.h".  Create new ASpect structure.  
#	16-Aug-94 by RAShaw	Make use of JDE's "aspect" memory structure 
#	28-Aug-94 by RAShaw	Make use of new database memory structure 
#	30-Aug-94 by RAShaw	Alter DB_READ calling sequence to include 
#				offset in running component number.   
#	24-Oct-94 by JDE	Change component list semantics & ASpect struct 
#	26-Apr-95 by RAShaw	Call as_clinit to fetch CL parameters; this 
#				used to be done in as_alloc(). 
#	27-Apr-95 by RAShaw	Use pointer to component list instead.
#	28-Apr-95 by RAShaw	Use component list managment routines.
#	24-May-95 by RAShaw	Change task name.
#	 5-Jun-95 by RAShaw	Support dispersion/flux units for output DB.
#	23-Jun-95 by RAShaw	Support input dispersion domains from a list.
#	10-Jul-95 by RAShaw	Check for missing DB components; call do_fit w/
#				size of component list rather than no_components
#				given in last DB header.
#	14-Jul-95 by RAShaw	Call do_fit w/component/spectrum list pointers.
#	27-Jul-95 by RAShaw	Loop over list of fit domains.
#	11-Aug-95 by RAShaw	Revise calling sequence to spmap.
#	12-Oct-95 by RAShaw	Fixed mem allocation bug for domain array. 
#	 5-Nov-95 by RAShaw	Accomodate new FIT structure.
#	 7-Nov-95 by RAShaw	Accomodate multiple input DB tables.
#	13-Dec-95 by RAShaw	Add explicit conversion of CP units; change 
#				calling args of db_write to use compon. list.
#	15-Dec-95 by RAShaw	Use new utility function to fetch desired 
#				units from the CL. 

procedure t_autospec()

#  Local variables:
pointer	as			# descriptor for ASpect memory structure
pointer	db			# DataBase descriptor
char	dbtab[SZ_FNAME]		# Name of component input table
pointer	domain_arr		# array of fit domain upper/lower limits
char	domain_list[SZ_FNAME]	# list of dispersion domains for fit
char	domain_stg[SZ_LINE]	# string of dispersion domains for fit
int	dunit, funit		# dispersion/flux unit enumerated type
double	dx1, dx2		# generic
#char	errmsg[SZ_LINE]		# text of error messages
pointer	ft			# fit object
int	i			# loop index
#char	input[SZ_FNAME]		# Name of spectrum input file
pointer	list			# list from input filename template
int	nft_reg			# number of fit intervals
bool	new			# write new output table?
pointer	o			# domain object
char	odbtab[SZ_FNAME]	# Name of component output table
#pointer	sp			# single spectrum object
int	status			# return status
char	input_list[SZ_LINE]	# string of input file names

#  Functions used:
pointer	as_alloc()		# initialize ASpect memory structure
extern	chispec()		# chi-square merit function
pointer	db_alloc()		# open DB table & initialize memory structure
int	do_fit()		# execute the fit, return status
int	get_domains()		# read dispersion domains from a string
int	imtgetim()		# get next image name from template
pointer	imtopenp()		# Open CL parameter list
int	next_domain()		# fetch next domain from list  
pointer	set_domain()		# allocate domain list parameter 
pointer	spmap()			# map a spectrum descriptor
int	tbtacc()		# test for existance of table

extern	fit_report()		# Report on progress of fitting.

errchk	as_alloc, as_free, clg_units, db_alloc, db_free, db_read, db_write
errchk	do_fit, set_domain, spl_add, spmap

begin

	# Set output DB flux/dispersion units.  If these are invalid, 
	# there's no point in continuing.
	call clg_units (dunit, funit)

	## Allocating ASpect memory structure
	as = as_alloc()
	call as_clinit(as)

	## Fetch input spectrum file names
	##list = imtopenp ("input")
        call clgstr ("input", input_list, SZ_LINE)
	call pr_sets (SP_INPARS(as), "flux", input_list)
	SP_LIST(as) = spmap (SP_INPARS(as))

	#  Read input spectra from disk
#	Note: the following two lines don't work, as CL list structured 
#		parameters are evidently not supported.  
#	call pr_sets (SP_INPARS(as), "flux", list)
#	SP_LIST(as) = spmap (SP_INPARS(as))

	#while (imtgetim (list, input, SZ_FNAME) != EOF) {
#	    call pr_sets (SP_INPARS(as), "flux", input)
#	    iferr ( sp = spmap (SP_INPARS(as)) ) {
#		call sprintf (errmsg, SZ_LINE, "Unable to open spectrum: %s")
#		    call pargstr (input)
#		call erract (EA_WARN)
#		call error (FIL_ACCESS, errmsg)
#	    }
#
#	    #  Merge the input spectrum list with the main list.
#	    do i = 1, SPL_N(sp) 
#		call spl_add (SP_LIST(as), SPL_A(sp,i), true)
#
#	    call spl_free (sp, false)
#	}

	#  Convert the input list of spectra to internal units.
	call spl_cv_units (SP_LIST(as), PHOTLAM, ANGSTROM)

	#  Fetch shortest/longest wavelengths.
	call spl_wave_domain (SP_LIST(as), dx1, dx2)
	call as_add_range (as, dx1, dx2)

	#  Examine spectrum objects
	if (DEBUG) {
	    do i = 1, SPL_N(SP_LIST(as))
	    	call sp_examine ( SPL_A(SP_LIST(as),i) )
	}

	## Now load components from input DB tables

	#  Open input DB table list and fetch components from each DB table
	list = imtopenp ("db_intab")
	while (imtgetim (list, dbtab, SZ_FNAME) != EOF) {
	    call get_cl_from_db (CMP_LIST(as), dbtab)
	}

	# Initialize user-function attributes from the ASpec object defaults.
	# call cp_uf_init (as, CMP_LIST(as))

	# Convert units of components to internal
	call cpl_cv_units (CMP_LIST(as), PHOTLAM, ANGSTROM)

	# Examine component objects
	if (DEBUG) {
	    do i = 1, CPL_N(CMP_LIST(as))
	    	call cp_examine ( CPL_A(CMP_LIST(as),i) )
	}

	## Enter loop to perform fits over requested domains

	call calloc (domain_arr, MAX_DOMAINS*2, TY_DOUBLE)

	#  Fetch fit domains from the CL
	call clgstr ("fit_domains", domain_list, SZ_FNAME)
	o = set_domain (domain_list)
	while (next_domain (o, domain_stg, SZ_LINE) != EOF) {
	    call printf ("Fit domain is: %s\n")
		call pargstr (domain_stg)

	    # Construct the fit control structure
	    ft = AS_FIT(as)
	    nft_reg = get_domains (domain_stg, Memd[domain_arr], WLOW(as), 
				   WHIGH(as), MAX_DOMAINS)
	    if (nft_reg <= 0) 
	    	call error (FIT_RANGE, "No valid fit domains")

	    # Initialize the FIT object with type, no. iterations, and 
	    # tolerances that reside in the global ASpect Fit object.  
	    call fit_init ( ft, SP_LIST(as), CMP_LIST(as), domain_arr, 
		nft_reg, FT_ALGORITHM(AS_FIT(as)), FT_ITER_LIMIT(AS_FIT(as)), 
		FT_DTOLER(AS_FIT(as)) )

	    # Execute the fit.
	    status = do_fit (ft, fit_report)

	    # Set the error calculation method and compute parameter errors.
	    call fit_seti (ft, FP_ERR_METHOD, CURV_MATRIX)
	    call do_errors (ft, chispec)

	    # Report error status, if any
	    if (status == ITER_EXCEEDED)
	    	call printf ("WARNING: Max. iterations exceeded during fit\n")
	    call printf ("\n")
	}

	call free_domain (o)
	call mfree (domain_arr, TY_DOUBLE)

	## Now write components to the output DB table.

	#  Open output DB table 
	call clgstr ("db_outtab", odbtab, SZ_FNAME)
	if (tbtacc (odbtab) == YES) {
	    new = false
	    db = db_alloc (odbtab, READ_WRITE)

	} else {
	    new = true
	    db = db_alloc (odbtab, NEW_FILE)
	    call db_set_units (db, funit, dunit)
	}

	call db_write (db, CMP_LIST(as), new) 

	## Close DB table, release memory & quit
	if (db != NULL)
	    call db_free (db)
	call as_free (as)

end


#-------------------------------------------------------------------------------
#  GET_CL_FROM_DB --	Build component list (DB_CL) from disk DB file and 
#			append to input CL.  
#
#  Revision history:
#	14-Dec-95 by RAShaw	Initial SPP implementation

procedure get_cl_from_db (cl, db_table)

#  Calling arguments:
pointer	cl			#I/O: component list
char	db_table[ARB]		#I: name of input DB table of components 

#  Local variables:
pointer	db			# input DataBase descriptor
pointer	db_cl			# component list read from DB
char	errmsg[SZ_LINE]		# text of error messages
int	ncp_expected		# no. components expected in DB
int	ncp_db			# no. components actually read from DB

#  Functions used:
pointer	db_alloc()		# open DB table & initialize memory structure
int	db_get_ncmp()		# get no. components in a DB object
int	db_read()		# read components from DB table 

errchk	db_alloc, db_read

begin
	db_cl = NULL
	db = db_alloc (db_table, READ_WRITE)

	# Ensure no. components matches that expected. 
	ncp_expected = db_get_ncmp (db)
	if (ncp_expected <= 0) {
	    call sprintf (errmsg, SZ_LINE, "DataBase table %s has zero entries")
		call pargstr (db_table)
	    call error (DB_COMPON, errmsg)
	}

	iferr ( ncp_db = db_read (db, db_cl, 0) ) {
	    call sprintf (errmsg, SZ_LINE, "Error reading DataBase table: %s")
		call pargstr (db_table)
	    call erract (EA_WARN)
	    call error (DB_COMPON, errmsg)
	}
	if (ncp_expected != ncp_db) {
	    call sprintf (errmsg, SZ_LINE, "Error reading DataBase table: %s")
		call pargstr (db_table)
	    call error (DB_COMPON, errmsg)
	}

	# Merge the input component list with the main list.
	call cpl_append (db_cl, cl, true)
	call cpl_free (db_cl, false)

	if (db != NULL)
	    call db_free (db)
end


#-------------------------------------------------------------------------------
#  FIT_REPORT -- Report fit progress to STDOUT. 

procedure fit_report (iter, chisq, force_flush)

#  Calling arguments:
int	iter			#I: The current iteration count.
real	chisq			#I: The current chisquare value.
bool	force_flush		#I: Force flushing the output?

begin
	call printf ("   iterations = %d chisq = %g\r")
	    call pargi (iter)
	    call pargr (chisq)

#	if (force_flush) 
	    call flush (STDOUT)

end
#-------------------------------------------------------------------------------
# End of fit_report
#-------------------------------------------------------------------------------