path: root/noao/twodspec/multispec/fitgauss5.x

include <fset.h>
include "ms.h"

# FITGAUSS5 -- Procedures used in fitting the GAUSS5 model.
#
# G5_FIT1 -- Fitting algorithm #1.
# G5_FIT2 -- Fitting algorithm #2.
# G5_FIT -- Fit the selected parameters for the best RMS value.
# SET_VERBOSE -- Verbose output.

########################################################################
.helpsys g5fit1 Jul84 MULTISPEC
.ih
NAME
G5_FIT1 -- Fitting algorithm #1.
.ih
DESCRIPTION
This algorithm fits the selected parameters simultaneously.
The parameters are selected by the parameter array (which of the 5
model parameters in a profile are to be fit) and the spectra range
array defined in fitgauss5.com.  These two arrays are used to generate
the fitparms array with the routine set_fitparams.  The fitparams array
controls the parameters fit by G5_FIT.
.ih
PARAMETERS
The model parameter values which are part of the MULTISPEC data structure,
the data array, the model array, the model profiles, and the profile
ranges must be initialized.  Other parameters for this procedure are
input via the common in file fitgauss5.com.  These include the spectra
to be fit and the parameters array.
.ih
OUTPUT
The returned data are the final parameter values, the model and profiles
arrays and the Y/N function value indicating if the RMS fit has been improved.
.endhelp
########################################################################

int procedure g5_fit1 (ms, data, model, profiles, ranges, lower, len_profile)

pointer	ms				# MULTISPEC database data
real	data[ARB]			# Line of image pixels
real	model[ARB]			# Line of model pixels
real	profiles[ARB]			# Model profiles
real	ranges[ARB]			# Origins of model profiles
real	lower				# Profile origin
int	len_profile			# Length of a model profile

int	improved
real	rms
pointer	sp, fitparams

int	g5_fit()
real	armsrr()

include	"fitgauss5.com"

begin
	# Calculate the initial RMS.  The parameter values are only changed
	# if the new RMS is less than this value.
	rms = armsrr (data, model, MS_LEN(ms, 1))
	call g5_prnt3 (rms)

	# Allocate and set the fitparams array.
	call smark (sp)
	call salloc (fitparams, MS_NSPECTRA(ms) * MS_NGAUSS5, TY_REAL)
	call set_fitparams (spectra, parameters, MS_NSPECTRA(ms), MS_NGAUSS5,
	    Memr[fitparams])

	# Call the fitting program once to simultaneously minimize the RMS.
	improved = g5_fit (ms, data, model, profiles, ranges, Memr[fitparams],
	    lower, len_profile, rms)

	call sfree (sp)
	return (improved)
end

###############################################################################
.helpsys g5fit2 Jul84 MULTISPEC
.ih
NAME
G5_FIT2 -- Fitting algorithm #2.
.ih
DESCRIPTION
This algorithm begins by fitting the parameters I0, X0, and S0
simultaneously.  Note that the values of S1 and S2 are used but are
kept fixed.  Next the parameters S0 and S1 (the shape) are fit simultaneously
keeping I0, X0, and S2 fixed followed by fitting I0 and X0 while
keeping S0, S1, and S2 (the shape) fixed.  If either of these fits
fails to improve the RMS then the algorithm terminates.
Also, if after the two steps (the fit of S0 and S1 followed by the fit
of I0 and X0), the RMS of the fit has not improved by more than the
user specified factor the algorithm also terminates.
.ih
INPUT
The model parameter values which are part of the MULTISPEC data structure,
the data array, the model array, the model profiles, and the profile
ranges must be initialized.  Other parameters for this procedure are
input via the common in file fitgauss5.com.  These include the spectra
to be fit, the parameters array (used as a working array), and the RMS
stopping factor.
.ih
OUTPUT
The returned data are the final parameter values, the model and profiles
arrays and the Y/N function value indicating if the RMS fit has been improved.
.endhelp
##############################################################################

int procedure g5_fit2 (ms, data, model, profiles, ranges, lower, len_profile)

pointer	ms				# MULTISPEC database data
real	data[ARB]			# Line of image pixels
real	model[ARB]			# Line of model pixels
real	profiles[ARB]			# Model profiles
real	ranges[ARB]			# Origins of model profiles
real	lower				# Profile origin
int	len_profile			# Length of a model profile

int	improved, fit
real	rms, rms_old
pointer	sp, fitparams

int	g5_fit()
real	armsrr()

include	"fitgauss5.com"

begin
	# Calculate the initial RMS.  The parameter values are only changed
	# if the new RMS is less than this value.
	rms = armsrr (data, model, MS_LEN(ms, 1))
	call g5_prnt3 (rms)

	# Allocate the fitparams array.
	call smark (sp)
	call salloc (fitparams, MS_NSPECTRA(ms) * MS_NGAUSS5, TY_REAL)

	# Fit the parameters I0, X0, and S0.
	parameters[I0_INDEX] = YES
	parameters[X0_INDEX] = YES
	parameters[S0_INDEX] = YES
	parameters[S1_INDEX] = NO
	parameters[S2_INDEX] = NO
	call set_fitparams (spectra, parameters, MS_NSPECTRA(ms),
	    MS_NGAUSS5, Memr[fitparams])

	# Call the fitting procedure to minimze the RMS.
	improved = g5_fit (ms, data, model, profiles, ranges,
	    Memr[fitparams], lower, len_profile, rms)

	# Two step fitting algorithm consisting of a fit to S0 and S1 followed
	# by a fit to I0 and X0.  This loop terminates when either one
	# of the fits fails to improve the RMS or the RMS has improved
	# by less than factor after the second step (the I0, X0 fit).
	repeat {
	    rms_old = rms

	    # Fit S0 and S1.
	    parameters[I0_INDEX] = NO
	    parameters[X0_INDEX] = NO
	    parameters[S0_INDEX] = YES
	    parameters[S1_INDEX] = YES
	    call set_fitparams (spectra, parameters, MS_NSPECTRA(ms),
		MS_NGAUSS5, Memr[fitparams])
	    fit = g5_fit (ms, data, model, profiles, ranges,
		Memr[fitparams], lower, len_profile, rms)
	    if (fit == NO)
		break
	    improved = YES

	    # Fit I0 and X0.
	    parameters[I0_INDEX] = YES
	    parameters[X0_INDEX] = YES
	    parameters[S0_INDEX] = NO
	    parameters[S1_INDEX] = NO
	    call set_fitparams (spectra, parameters, MS_NSPECTRA(ms),
		MS_NGAUSS5, Memr[fitparams])
	    fit = g5_fit (ms, data, model, profiles, ranges,
		Memr[fitparams], lower, len_profile, rms)
	    if (fit == NO)
	        break

	    if (rms > (1 - factor) * rms_old)
		break
	}

	call sfree (sp)
	return (improved)
end


##############################################################################
.helpsys g5fit Jul84 MULTISPEC
.ih
NAME
G5_FIT -- Basic parameter fitting procedure.
.ih
INPUT
The input data are the data array to be fit and the initial model
parameters (part of the MULTISPEC data structure), the model array
and model profiles (with the profile ranges array) corresponding to the
initial model parameters, and the RMS of the model relative to the data.
The parameters to be fit are selected by the fitparams array.
Parameters controlling the fitting process are input to this procedure
via the common block in the include file fitgauss5.com.  These parameters are
the RMS stopping factor and parameters controlling the smoothing of the
shape parameters.
.ih
OUTPUT
The returned data are the final parameter values, the model and profiles
arrays and the Y/N function value indicating if the RMS fit has been improved.
.ih
DESCRIPTION
The best RMS fit is obtained by iteration.  Correction vectors for the
parameters being fit are obtained by the simultaneous banded matrix
method in the procedure solve.  Heuristic constraints and smoothing
are applied to the solution and then the RMS of the new fit to the
data is calculated.  New parameter corrections are computed until the RMS of
the fit fails to improve by the specified factor.
.endhelp
############################################################################

int procedure g5_fit (ms, data, model, profiles, ranges, fitparams, lower,
    len_profile, rms)

pointer	ms				# MULTISPEC data structure
int	fitparams[ARB]			# Model parameters to be fit
real	data[ARB]			# Data line to be fit
real	model[ARB]			# Model line
real	profiles[ARB]			# Model profiles
real	ranges[ARB]			# Profile ranges
real	lower				# Lower limit of profiles
int	len_profile			# Length of profiles
real	rms				# RMS of fit

int	improved
int	len_line, nspectra, nparams
real	rms_next, rnorm
pointer	sp, last_i0, last_x0, last_s0, last_s1, last_s2
pointer solution, sol_i0, sol_x0, sol_s0, sol_s1, sol_s2

real	armsrr()

include	"fitgauss5.com"

begin
	# Set array lengths.
	len_line = MS_LEN(ms, 1)
	nspectra = MS_NSPECTRA(ms)
	nparams = MS_NGAUSS5

	# Allocate working memory to temporarily save the previous parameter
	# values and to hold the correction vector.
	call smark (sp)
	call salloc (last_i0, nspectra, TY_REAL)
	call salloc (last_x0, nspectra, TY_REAL)
	call salloc (last_s0, nspectra, TY_REAL)
	call salloc (last_s1, nspectra, TY_REAL)
	call salloc (last_s2, nspectra, TY_REAL)
	call salloc (solution, nspectra * nparams, TY_REAL)

	# Offsets in the solution array for the various parameters.
	sol_i0 = solution + (I0_INDEX - 1) * nspectra
	sol_x0 = solution + (X0_INDEX - 1) * nspectra
	sol_s0 = solution + (S0_INDEX - 1) * nspectra
	sol_s1 = solution + (S1_INDEX - 1) * nspectra
	sol_s2 = solution + (S2_INDEX - 1) * nspectra

	improved = NO
	repeat {
	    # Store the last parameter values so that if the parameter values
	    # determined in the next iteration yield a poorer RMS fit to
	    # the data the best fit parameter values can be recovered.

 	    call amovr (PARAMETER(ms,I0,1), Memr[last_i0], nspectra)
 	    call amovr (PARAMETER(ms,X0,1), Memr[last_x0], nspectra)
 	    call amovr (PARAMETER(ms,S0,1), Memr[last_s0], nspectra)
 	    call amovr (PARAMETER(ms,S1,1), Memr[last_s1], nspectra)
 	    call amovr (PARAMETER(ms,S2,1), Memr[last_s2], nspectra)

	    # Determine a correction solution vector for the selected
	    # parameters simultaneously, apply heuristic constraints to the
	    # solution vector, apply the correction vector to the parameter
	    # values, and smooth the shape parameters if requested.

	    # Find a least squares correction vector.
	    call solve (ms, data, model, fitparams, profiles, ranges,
		len_line, len_profile, nspectra, nparams, Memr[solution], rnorm)

	    # Apply constraints to the correction vector.
	    call constrain_gauss5 (ms, Memr[solution], nspectra, nparams)

	    # Add the correction vector to the parameter vector.
	    call aaddr (PARAMETER(ms,I0,1), Memr[sol_i0], PARAMETER(ms,I0,1),
		nspectra)
	    call aaddr (PARAMETER(ms,X0,1), Memr[sol_x0], PARAMETER(ms,X0,1),
		nspectra)
	    call aaddr (PARAMETER(ms,S0,1), Memr[sol_s0], PARAMETER(ms,S0,1),
		nspectra)
	    call aaddr (PARAMETER(ms,S1,1), Memr[sol_s1], PARAMETER(ms,S1,1),
		nspectra)
	    call aaddr (PARAMETER(ms,S2,1), Memr[sol_s2], PARAMETER(ms,S2,1),
		nspectra)

	    # Smooth the shape parameters.
	    if (smooth[S0_INDEX] == YES)
		call ms_smooth (PARAMETER(ms, X0, 1), PARAMETER(ms, S0, 1))
	    if (smooth[S1_INDEX] == YES)
		call ms_smooth (PARAMETER(ms, X0, 1), PARAMETER(ms, S1, 1))
	    if (smooth[S2_INDEX] == YES)
		call ms_smooth (PARAMETER(ms, X0, 1), PARAMETER(ms, S2, 1))

	    # Calculate new model profiles and new model data line.
	    # Determine the RMS fit of the new model to the data.
	    # If the change in the RMS is less than factor times the
	    # previous RMS the interation is terminated else the improvement
	    # in the RMS is recorded and the next iteration is begun.

	    # Set new model profiles.
	    call mod_gauss5 (ms, lower, profiles, ranges, len_profile, nspectra)

	    # Set new model line from the profiles.
	    call set_model (ms, model, profiles, ranges, len_line,
		len_profile, nspectra)

	    # Calculate the RMS of the new model.
	    rms_next = armsrr (data, model, len_line)

	    # Check to see if the RMS is improved enough to continue iteration.
	    if ((rms - rms_next) < factor * rms) {

		# The RMS has not improved enough to continue iteration.

		if (rms_next < rms) {
		    # Keep the latest parameter values, profiles, and model
		    # because the new RMS is lower than the previous RMS.
		    # Record the improvement.
		    rms = rms_next
		    improved = YES
		    call g5_prnt3 (rms)

		} else {
		    # Restore the parameter values, profiles, and model to 
		    # previous values because the new RMS is higher.
 	    	    call amovr (Memr[last_i0], PARAMETER(ms,I0,1), nspectra)
 	    	    call amovr (Memr[last_x0], PARAMETER(ms,X0,1), nspectra)
 	    	    call amovr (Memr[last_s0], PARAMETER(ms,S0,1), nspectra)
 	    	    call amovr (Memr[last_s1], PARAMETER(ms,S1,1), nspectra)
 	    	    call amovr (Memr[last_s2], PARAMETER(ms,S2,1), nspectra)
	    	    call mod_gauss5 (ms, lower, profiles, ranges, len_profile,
			nspectra)
	    	    call set_model (ms, model, profiles, ranges, len_line,
			len_profile, nspectra)
		}

		# Exit the iteration loop.
		break

	    } else {

		# The RMS has improved significantly.  Record the improvement
		# and continue the iteration loop.

	        rms = rms_next
	        improved = YES
		call g5_prnt3 (rms)
	    }
	}

	call sfree (sp)
	return (improved)
end


# G5_SET_VERBOSE -- Output procedures for verbose mode.

procedure g5_set_verbose (verbose)

bool	verbose
bool	flag

# entry g5_prnt1 (image, naverage, track, start)
char	image[1]
int	naverage
bool	track
int	start

# entry g5_prnt2 (line, data, len_data)
int	line, len_data
real	data[1]
real	rms, data_rms

real	armsr()
include	"fitgauss5.com"

begin
	# Toggle verbose output.
	flag = verbose
	if (flag)
	    call fseti (STDOUT, F_FLUSHNL, YES)
	else
	    call fseti (STDOUT, F_FLUSHNL, NO)
	return

entry g5_prnt1 (image, naverage, track, start)

	# Print the values of the various task parameters.

	if (!flag)
	    return

	call printf ("\nMULTISPEC Model Fitting Program\n\n")
	call printf ("Image file being modeled is %s.\n")
	    call pargstr (image)
	call printf ("Average %d lines of the image.\n")
	    call pargi (naverage)
	call printf ("Fitting algorithm %d.\n")
	    call pargi (algorithm)
	if (algorithm == 1) {
	    if (parameters[I0_INDEX] == YES)
	        call printf ("Fit intensity scales.\n")
	    if (parameters[X0_INDEX] == YES)
	        call printf ("Fit spectra positions.\n")
	    if (parameters[S0_INDEX] == YES)
	        call printf ("Fit spectra widths.\n")
	    if (parameters[S1_INDEX] == YES)
	        call printf ("Fit model parameter s1.\n")
	    if (parameters[S2_INDEX] == YES)
	        call printf ("Fit model parameter s2.\n")
	}
	if (track) {
	    call printf ("Track model from line %d.\n")
		call pargi (start)
	}
	call printf (
	    "Iterate model until the fit RMS decreases by less than %g %%.\n\n")
	    call pargr (factor * 100)

	return

entry g5_prnt2 (line, data, len_data)

	# Print the image line being fit and the data RMS.
	if (flag) {
	    call printf ("Fit line %d:\n")
	        call pargi (line)
	    data_rms = armsr (data, len_data)
	    call printf ("  Data RMS = %g\n")
	        call pargr (data_rms)
	}
	return

entry g5_prnt3 (rms)

	# Print the RMS of the fit and the ratio to the data RMS.
	if (flag) {
	    call printf ("  Fit RMS = %g  Fit RMS / Data RMS = %g\n")
	        call pargr (rms)
	        call pargr (rms / data_rms)
	}
end