; Copyright (c) 1998-2005 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_SVD_spectra
;
;LAST CHANGED: ----------------------------------- 03-feb-05
;
;PURPOSE:
;	This procedure uses SVD to fit a spectrum to a set of reference spectra.
; AX_SVD_spectra will run stand alone (prompting user for AXIS format files) or from AXIS.
;
;CATEGORY:
;	STAND ALONE: spectral analysis
;
;CALLING SEQUENCE:
;	AX_SVD_spectra[SPC, axis=axis, help = help, coeff = coeff, set_start = set_start, pos_only=pos_only ]
;
;CALLED FROM AXIS:
;	->{Spectra}->Curve Fit~SVD_spectra

;INPUTS:
;	SPC		spectrum to be fit in axis 1d structure form
;
;KEYWORDS:
; AXIS 		called from AXIS widget
; COEFF		filename with spectral information
; HELP		print help text
; SET_START	allow users to set starting values
; POS_ONLY	force only positive coefficients
; LIMITS	[min, max] 2-vector defining lower and upper limits
;
;OUTPUTS: 	No explicit outputs.
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;PROCEDURE:
;	The reference spectra are read in and the coefficient matrix generated.
; SVD is then called to get the fittin coefficients
;
;MODIFICATION HISTORY:
; (03-feb-04 aph) first version - adapted from ax_curvfit (20-oct-04 version )
;-

pro ax_SVD_spectra, spc, axis=axis, help=help, coeff = coeff, set_start=set_start, pos_only=pos_only

@axis_com

IF keyword_set(help) THEN BEGIN
    print,'AX_SVD_spectra'
    print,'Fits a spectrum to a set of reference spectra using SVD'
    print, 'Uses AXIS format spectra files (*.txt) as input/output'
    print, ' KEYWORDS: '
    print, '	AXIS   = if on, then called from AXIS widget'
  	print, '	HELP   = print this how-to-use statement'
  	print, '	SET_START = allow user to define starting parameters'
  	print, '	POS_ONLY = force fit coefficients to be positive'
  	print, '    LIMITS = define bounds for fit coefficients'
  	print, '    SWAP = if on, make too low values upper and vice-versa'
    return
ENDIF

; assume AXIS is running (therefore may have called ax_curvfit)
; either when AXIS keyword is supplied or if any widget active
if  keyword_set(axis) then axis_on = 1 else axis_on = widget_info(/active)
print, ' '
print, ' Curve fit analysis by CG_Optimize'
print, ' -------------------------------------------'

;--------- get the spectrum to be fit
if n_tags(spc) EQ 0 then begin
	text =  ' Select spectrum to fit'
	print, text
	widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
	spc = spc_load()
	if n_tags(spc) EQ 0 then return		; bail if user says cancel
ENDIF


; -------------- define coefficents at energies of images -------------
IF NOT keyword_set(coeff) THEN BEGIN
	t = ' '
    t = dialog_message('Read fit parameter file ?', /question)
	if t(0) EQ 'Yes' then begin
	    par_file = pickfile2(title = ' Select fit parameter file', $
	           filter='*.par', /LPATH, DEFPATH=defpath )
	    if strlen(par_file(0)) NE 0 then coeff = par_file
	endif
ENDIF
IF keyword_set(coeff) THEN BEGIN
; ------- read spectral file info from a file
	check = findfile(coeff)
	if strlen(check(0)) NE 0 then pars = ax_par_load(coeff) else goto, get_user_info
	if n_tags(pars) EQ 0 then begin
		goto, get_user_info
	endif else begin
		ncomp = pars.n
		comp_names = pars.names
		comp_files = pars.files
	endelse

; ------ if fails then ask user to locate reference files one-by-one
ENDIF ELSE BEGIN
 get_user_info:
	if axis_on then ncomp = get_num(prompt = 'Number of components',val=ncomp, group=axis_ID) $
	   else ncomp = get_num(prompt = 'Number of components',val=ncomp)
	if ncomp LE 0 then return
	comp_names = strarr(ncomp)
	comp_files = strarr(ncomp)
; ------- read coefficent information from spectra files
	for i = 0, ncomp-1 do begin
	 	text = 'Spectral file for component ' + strcompress(string(fix(i+1)))
		if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
		comp_files(i) = pickfile2(title = ' Spectrum of component '+ strtrim(string(i),2), $
	           filter='*.txt', /LPATH, DEFPATH=defpath )
	    if comp_files(i) EQ '' then return
		tmp = spc_load(file = comp_files(i))
;; ------ define name - keep short as used in map name
		text = 'short name of component ' + strcompress(string(fix(i+1)))
		if axis_on then comp_names(i) = get_text(prompt = text, val = tmp.dl, group = axis_ID) $
		    else comp_names(i) = get_text(prompt = text, val = tmp.dl)
		comp_names(i) = strtrim(string(comp_names(i)),2)
	endfor

; ------ save coefficient array for later use (optional)
	par_file = pickfile2(filter='*.par', /write, path = DefPath, $
	           title = 'Name of fit parameter file')
	if strlen(par_file(0)) NE 0 then ax_par_save,par_file,ncomp,comp_names,comp_files
ENDELSE

; ------------- check for out-of-order or non-monotonic data sets -----------------------------
; (aph) 30-dec-03 - if there are two equal x-values in a spectrum scale this leads to NaN fits due to
;                   failure if the IDL INTERPOL algorithm (NB requirement of monotonic data is noted in the
;                   header for INTERPOL.PRO (IDL5.2) and in the HELP file

 spc = ax_sort_mono(spc,/axis)

; read in files and interpolate to the energy scale of the SPECTRUM
A=fltarr(ncomp, n_elements(spc.x))
CurBuf = 0
for i = 0, ncomp-1 do begin
	tmp = spc_load(file=comp_files(i))
		tmp = ax_sort_mono(tmp,/axis)			; FORCE MONOTONIC !!
	; ----- interpolate to the same energies
	a(i,*) = interpol(tmp.d, tmp.x, spc.x)		; FORCE INTERPOLATION !!

; ----------- plot reference spectra if axis on
	if axis_on then begin
		CurBuf = i+1
		stmp ={t:'1d', x: spc.x, d: a(i,*), dn: a(i,*), dl: comp_names(i)}
	   	HANDLE_value, Data(CurBuf),stmp,/set
		Label(CurBuf) = stmp.dl
		Plotbuf, CurBuf
	endif
endfor
; -------------- define starting coefficents of fit -------------
; let USER CHOOSE start value of unknown mix coeff
	xp = fltarr(ncomp)
	if keyword_set(set_start) then begin
		if axis_on THEN BEGIN
			for i = 0, ncomp-1 do $
				xp(i) = get_num(prompt=comp_names(i) + ' start value', val=0, group = axis_ID)
		endif else for i = 0, ncomp-1 do $
				xp(i) = get_num(prompt=comp_names(i) + ' start value', val=0)
	endif

; ------------ execute AX_CurvFit procedure --------------------

numit = 1000

; ---------- LET USER CHOOSE tolerance
tolerance = MACHAR()
TOL=  0.01*SQRT(tolerance.eps)
Tol = 1.e-20
if axis_on THEN tol = get_num(prompt='tolerance', val=tol, group = axis_ID) else $
	tol = get_num(prompt='tolerance', val=tol)

; ---- optional force all positive coefficients, check if not set on call
IF NOT keyword_set(pos_only) then begin
	t = dialog_message('Force only positive values ?', /question, /default_no)
	if t(0) EQ 'Yes' then pos_only = 1 else pos_only = 0
ENDIF

; ---- optional limits
IF NOT keyword_set(limits) AND pos_only NE 1 then begin
	t = dialog_message('Set limits ?', /question, /default_no)
	if t(0) EQ 'Yes' then begin
		limits = fltarr(2)
		if axis_on THEN limits(0) = get_num(prompt='lower limit', group = axis_ID) else $
			limits(0) = get_num(prompt='lower limit')
		if axis_on THEN limits(1) = get_num(prompt='upperr limit', group = axis_ID) else $
			limits(1) = get_num(prompt='upper limit')
		if limits(0) GT limits(1) then begin		; force lower/upper
			tt = limits(0) & limits(0) = limits(1) & limits(1)= tt
		endif
		tt = dialog_message('Swap if hit limits ?', /question, /default_no)
		if tt(0) EQ 'Yes' then swap = 1 else swap = 0
	endif
ENDIF ELSE begin
	swap = 0 & limits = 0
ENDELSE

; carry out CG_Optimize ---------
widget_control,/hourglass
xp = CG_OPTIMIZE(xp,'cgex_func', 'cgex_dfunc', LINMIN='cgex_linmin', $
		num_iter=numit, tolerance = tol, pos_only = pos_only, $
		limits = limits, swap = swap, $
	    OBJECTIVE=obj, A=A, b=spc.d, _extra=e)

; compute fit and residuals
fit = spc
fit.d = TRANSPOSE(a##xp)
fit.dl = 'fit to data'
res = spc
res.d = spc.d - fit.d
res.dl = 'residual of fit'

print, 'CG_Optimize:  number of iterations = ', numit
t = moment(res.d, sdev = sdev)
print, ' Standard deviation = ', sdev
print, ' Fit components', xp

; ---- output spectral components
if axis_on then begin
; report the fit results
	text = string(xp)
	widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
; store fit components and residual in axis buffers
	offset = 1  & td = fltarr(n_elements(spc.d))
	for i = 0, ncomp-1 do begin
		widget_control,/hourglass
		CurBuf = i + offset
		s = spc
		s.d = xp(i)*a(i,*)
		td = td + s.d
		s.dl = strtrim(string(xp(i)),2) + ' * ' + comp_names(i)
	   	HANDLE_value, Data(CurBuf),s,/set
		Label(CurBuf) = s.dl
		Plotbuf, CurBuf
	endfor

; output fit spectrum
	CurBuf = CurBuf + 1
	s.d = td
	s.dl = strtrim(string(ncomp),2) + ' FIT of ' + spc.dl
	HANDLE_value, Data(CurBuf),s,/set
	Label(CurBuf) = s.dl
	Plotbuf, CurBuf

; place residuals in buffer 0
	CurBuf = 0
	HANDLE_value, Data(CurBuf),res,/set
	Label(CurBuf) = res.dl
	Plotbuf, CurBuf
endif else begin

; ---- plot fit if run in standalone mode
	ymin = min([spc.d,fit.d, res.d], max = ymax)
	splot, spc, yrange = [ymin, ymax]
	splot, fit, /o
	splot, res, /o
endelse

return

end

