; Copyright (c) 1998-2013 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_regress_spectra
;
;LAST CHANGED: ----------------------------------- 31-Dec-12
;
;PURPOSE:
;	This procedure uses linear regression to fit a spectrum to a set of reference spectra.
; AX_REGRESS_spectra will run stand alone (prompting user for AXIS format files) or from AXIS.
;
;CATEGORY:
;	STAND ALONE: spectral analysis
;
;CALLING SEQUENCE:
;	AX_REGRESS_SPECTRA[SPC, axis=axis, help = help]
;
;CALLED FROM AXIS:
;	->{Spectra}->Curve Fit~linear regression

;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
;
;OUTPUTS: 	No ecoefflicit outputs.
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;PROCEDURE:
;	The reference spectra are read in and the coefficient matrix generated.
; REGRESS is then called to get the fittin coefficients
; qauality of fit is evaluated by Ftest, Mcorrelation and chi-square
; for F-test - larger numbers are better
; for Mcorrelation - numbers closer to +1.0 are better; values close ot 0 are bad
; for chi-square - smaller numbers are better
;
;MODIFICATION HISTORY:
; (03-feb-05 aph) first version - adapted from ax_curvfit (20-oct-04 version)
; (19-feb-05 aph) correct final display
; (22-jun-05 aph) fix display of results (axis_log); switch to avoid 'regress' in IDL5.2
; (17-Jan-12 aph) add Mcorrelation evaluator of statistical precision
; (31-Dec-12 aph) add header to log to track data
;-

pro ax_regress_spectra, spc, axis=axis, help=help, coeff = coeff
@axis_com
on_error,2

IF keyword_set(help) THEN BEGIN
    print,'AX_regress_spectra'
    print,'Fits a spectrum to a set of reference spectra using IDL REGRESS function'
    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'
    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 linear regression'
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)

; --------- construct appropriate arrays ---------------

; 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+1 then CurBuf = 0 else CurBuf=CurBuf+1   ; 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


; ------------ execute linear regression procedure --------------------

widget_control,/hourglass
y_vals=spc.d
;  ------- no weighting (equal weighting of points
weights=make_Array(n_elements(y_vals), /float, value=1.0)
; ------- Gaussian weighting ------
; weights = 1.0/STDDEV(y_vals)^2.   ; THIS DOES NOT WORK. aph 22-jun-05
if !version.release GE 6.0 then begin
        coeff = Regress(A,y_vals,chisq=chisquare, const=const, correlation=correl, $
              ftest=ftest, measure_errors=weights, mcorrelation=mcorr, status=stat, yfit=fity)
endif else begin
		axis_log, 'WARNING: Linear regression in IDL 5.2 often does not work'
        coeff = Regress(A,y_vals, weights,fity, const, /relative_weight, $
                  sigma,ftest,rdum, rmul, chisquare, stat)
endelse
if stat NE 0 then begin
	if stat EQ 1 then text = 'singular array - inversion is invalid'
	if stat EQ 2 then text= 'warning that a small pivot element was used. Significant accuracy was probably lost'
	axis_log, text
endif


;-------- compute residuals ----------------------
fit = spc
fit.d = fity
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 = '--------------------------------------'
axis_log, text
text = 'Curve fit by linear regress (stack fit):  HEADER = ' + spc.dl
axis_log, text
text = '--------------------------------------'
axis_log, text


; if !version.release GE 6.0 then begin
	text = ' F-test = '+ set_digits(ftest,3)
	axis_log, text
		text = ' M_correlation = '+ set_digits(mcorr,3)
	axis_log, text
	text = ' chi square = ' + set_digits(chisquare,3)
	axis_log, text
; endif

axis_log, ' Fit components'
for i = 0,ncomp-1 do begin
	text =  comp_names(i) + '   ' + set_digits(coeff(i),3)
	if axis_on then axis_log, text
endfor
text =  'const   ' + set_digits(const,3)
axis_log, text


; ---- output spectral components
if axis_on then begin
; store fit components and residual in axis buffers
	if SpectrumBuf GT ncomp+1 then CurBuf = 1 else CurBuf=SpectrumBuf+1   ; either plot early or late
	for i = 0, ncomp-1 do begin
		widget_control,/hourglass
		s = spc
		s.d = coeff(i)*a(i,*)
		y_min = min([y_min,s.d])
		y_max = max([y_max,s.d])
		s.dl = strtrim(string(coeff(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
; generate constant component
	s = spc
	s.d(0:n_elements(s.d)-1) = const
	y_min = min([y_min,s.d])
	y_max = max([y_max,s.d])
	s.dl = string(const)+'  constant'
   	HANDLE_value, Data(CurBuf),s,/set
	Label(CurBuf) = s.dl
	Plotbuf, CurBuf

; output fit spectrum
	CurBuf = CurBuf + 1
	if CurBuf EQ 10 then CurBuf=1
	y_min = min([y_min,fit.d])
	y_max = max([y_max,fit.d])
	fit.dl = 'FIT of ' + spc.dl
	HANDLE_value, Data(CurBuf),fit,/set
	Label(CurBuf) = fit.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+1 then CurBuf = 1 else CurBuf=CurBuf+1
for i =0,ncomp+1 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
