; Copyright (c) 1998-2013  A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_CURVFIT
;
;LAST CHANGED: ----------------------------------- 31-Dec-12
;
;PURPOSE:
;	This procedure executes the CG_Optimize procedure (Billy Loo @2000)
; to fit a spectrum to a set of reference spectra.
; AX_CURVFIT will run stand alone (prompting user for AXIS format files) or from AXIS.
;
;CATEGORY:
;	STAND ALONE: spectral analysis
;
;CALLING SEQUENCE:
;	AX_CURVFIT[SPC, axis=axis, help = help, coeff = coeff, set_start = set_start, pos_only=pos_only ]
;
;CALLED FROM AXIS:
;	->{Spectra}->Curve fFit

;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.
; CG_OPTIMIZE is then called to get the mixing coefficients
;
;MODIFICATION HISTORY:
; (13-feb-01 aph) first version
; (27-mar-01 aph) changed title of file input dialog
; (26-jul-01 aph) neatened code ; get pos_only working; explicitly include A, B in call
; (04-jun-03 aph) adapt for reading ref spectra names etc from ax_par_load & save to ax_par_save
; (30-dec-03 aph) force monotonic spectral data (source of NaN fits !!)
; (20-oct-04 aph) set defaults to <No>
; (03-feb-05 aph) add auto-plot of all fit at end
; (19-feb-05 aph) correct display at end
; (22-jun=05 aph) correct syntax for axis_log calls
; (13-dec-12 aph) display the header of the data with log

pro ax_CurvFit, 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_CurvFit'
    print,'Executes CG_Optimize procedure of Billy Loo to fit a spectrum'
    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))
SpectrumBuf=CurBuf   ; retain value of spectrum buffer
if CurBuf GT ncomp then CurBuf = 0   ; either plot early or late
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 = CurBuf +1
		if CurBuf EQ 10 then CurBuf=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'

y_min = min(res.d, max=y_max)

; ----------- report results --------------------

text = 'Curve fit:  HEADER = ' + spc.dl
axis_log, text
text = 'CG_Optimize: ' + string(fix(numit))+' loops'
axis_log, text
t = moment(res.d, sdev = sdev)
text = ' Std. dev. = ' + set_digits(sdev,3)
axis_log, text
print, ' Fit components'
for i = 0,ncomp-1 do begin
	text =  comp_names(i) + '  ' + set_digits(xp(i),3)  ;strtrim(string(xp(i)),2)
	axis_log, text
endfor
;text =  string('const   ', const)
;if axis_on then axis_log, text else print, text


; ---- output spectral components
if axis_on then begin
; store fit components and residual in axis buffers
	if SpectrumBuf GT ncomp then CurBuf = 1 else CurBuf=SpectrumBuf+1   ; either plot early or late
	td = fltarr(n_elements(spc.d))
	for i = 0, ncomp-1 do begin
		widget_control,/hourglass
		s = spc
		s.d = xp(i)*a(i,*)
		td = td + s.d
		y_min = min([y_min,s.d])
		y_max = max([y_max,s.d])
		s.dl = strtrim(string(xp(i)),2) + ' * ' + comp_names(i)
	   	HANDLE_value, Data(CurBuf),s,/set
		Label(CurBuf) = s.dl
		Plotbuf, CurBuf
		CurBuf = CurBuf+1
		if CurBuf EQ 10 then CurBuf=1
	endfor

; output fit spectrum
	s.d = td
	y_min = min([y_min,s.d])
	y_max = max([y_max,s.d])
	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

; --------- plot complete fit
CurBuf = SpectrumBuf   ; recall spectrum buffer
y_offset=0.05*(y_max-y_min)
Yrng=[y_min-y_offset,y_max+y_offset]
SetGraf,'MainImg'
HANDLE_VALUE, Data(CurBuf), tmp
	splot, tmp, color=bcolor(CurBuf), charsize=0.7, yrange = Yrng, $
		   thick=4, xstyle=1, ystyle=1, psym=Line_sym

if CurBuf GT ncomp then CurBuf = 1 else CurBuf=CurBuf+1
for i =0,ncomp do begin
	HANDLE_VALUE, Data(CurBuf), tmp
	splot, tmp, color=bcolor(CurBuf), /o, charsize=0.7, yrange = Yrng, $
		   thick=2, xstyle=1, ystyle=1, psym=Line_sym
	CurBuf = CurBuf+1
	if CurBuf EQ 10 then CurBuf=1
endfor
CurBuf=0
HANDLE_VALUE, Data(CurBuf), tmp
splot, tmp, color=bcolor(CurBuf), /o, charsize=0.7, yrange = Yrng, $
	   thick=2, xstyle=1, ystyle=1, psym=Line_sym

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

