; Modifications:
;  23-Jun-2004. Kim. Call mcurvefit with p_noeffect.  If fail_type is 3 then set those
;    params to fixed, and try again.  Save free_mask used in info param mcurvefit_free_mask.
;  25-Aug-2005, Kim.  Added /allow_allfixed in call to mcurvefit
;  30-Sep-2005, Kim. Added class='fit_function' in call to get(/fit_function) for speed,
;    and use class='fit_comp_manager' in call to get fit_comp_params for speed
;  26-Apr-2006, Kim. Fixed small bug with setting fixed params with p_noeffect
;  8-Aug-2006, Kim.  Previously two sources - fitrange and fit_function as src_index=1.
;   Now only source is fit_function.  Spex_fit_obj is stored as property so we can get what
;   we need from it, without having it be part of this source chain.
;   Also, in process_hook, pass in loop_index for outer loop over fitting, and only
;   save #iterations when loop_index is 0.
;  27-Feb-2008, Kim.  Supply format for printing fit results params, and added printing sigmas
;  7-May-2008, Kim. Initialize sigmas to 0. in process_hook in case no fitting is done in mcurvefit
;  19-Feb-2009, Kim. In process_hook, increased precision for printing chi2
;  3-Dec-2009, Kim.  Added saving mcurvefit_corr correlation matrix
;  20-Sep-2010, Kim.  Print full chi as well as reduced after fitting.  Also added quiet keyword to 
;    process_hook, and set show_progress based on quiet value.
;  27-Apr-2012, Kim.  Added saving mcurvefit_covar covariance matrix
;  01-Aug-2014, Kim. Added number of degrees of freedom to printed summary of fit
;
;------------------------------------------------------------------------------
function spex_fitalg_mcurvefit::init, $
	source=source, $
	_extra=_extra

control = spex_fitalg_mcurvefit_control()
info = { spex_fitalg_mcurvefit_info }

if not obj_valid( source ) then source = obj_new('fit_function')

ret = self->framework::init( control=control, $
				info=info, $
				source=source, $
        		_extra=_extra  )


;each fit alg will have its own version of a fit_function object
;self -> framework::set, source = obj_new('fit_function'), src_index=1

if keyword_set(_extra) then self -> set, _extra=_extra

return, ret

end

;------------------------------------------------------------------------------

;pro spex_fitalg_mcurvefit::set, this_class_only=this_class_only, _extra=_extra
;
;self -> framework::set, _extra=_extra, this_class_only=0

;end

;------------------------------------------------------------------------------

function spex_fitalg_mcurvefit::getalg, chisq=chisq, niter=niter, maxiter=maxiter, stop_msg=stop_msg, free_mask=free_mask
if keyword_set(chisq) then return, self -> get(/mcurvefit_chi2)
if keyword_set(niter) then return, self -> get(/mcurvefit_iter)
if keyword_set(stop_msg) then return, self -> get(/mcurvefit_fail_msg)
if keyword_set(maxiter) then return, self -> get(/mcurvefit_itmax)
if keyword_set(free_mask) then return, self -> get(/mcurvefit_free_mask)
return, -1
end

;------------------------------------------------------------------------------
pro spex_fitalg_mcurvefit::process_hook, xdata, ydata, errors, $
	result=result, $
	params=params, $
	sigmas=sigmas, $
	loop_index=loop_index, $
	quiet=quiet
	
checkvar, quiet, 0	
show_progress = ~quiet

; it's faster to do this the second way
;params = self -> getfit(/fit_comp_params)
;minarr = self -> getfit(/fit_comp_minima)
;maxarr = self -> getfit(/fit_comp_maxima)
;freearr = self -> getfit(/fit_comp_free_mask)

params = self -> get(/fit_comp_params, class='fit_comp_manager')
minarr = self -> get(/fit_comp_minima, class='fit_comp_manager')
maxarr = self -> get(/fit_comp_maxima, class='fit_comp_manager')
freearr = self -> get(/fit_comp_free_mask, class='fit_comp_manager')

; set min and max of non-free parameters equal to parameter value
fixix = where (freearr eq 0, nfixed)

rep_count = 0
while rep_count lt 2 do begin
	rep_count = rep_count + 1
	if nfixed gt 0 then begin
		minarr[fixix] = params[fixix]
		maxarr[fixix] = params[fixix]
	endif
  
  sigmas = params * 0.
  
	result = mcurvefit( xdata, $
		ydata, $
		errors, $
		params, $
		sigmas, $
		corr=corr, $
		covar=covar, $
		function_name='spex_fitalg_mcurvefit_eval', $
		itmax=self->get(/mcurvefit_itmax), $
		iter=iter, $
		tol=self->get(/mcurvefit_tol), $
		chi2=chi2, $
		/noderivative, $
		minarr=minarr, $
		maxarr=maxarr, $
		fail_type=fail_type, $
		quiet=quiet, $
		show_progress=show_progress, $
		p_noeffect=p_noeffect, $
		/allow_allfixed, $
		private=self )

	if (fail_type ne 3) or (rep_count gt 1) or (p_noeffect[0] eq -1) then break else begin
		fixix = nfixed gt 0 ? [fixix, p_noeffect] : p_noeffect
		nfixed = n_elements(fixix)
		freearr = freearr*0 + 1
		freearr[fixix] = 0
		message,'Trying fit again with parameters ' + arr2str(trim(p_noeffect)) + ' fixed, not free.', /cont
	endelse
endwhile

; mcurvefit calculates chisq divided by (# data points - # params).  We want chisq divided by
; (# data points - # free params -1 )
;npts = n_elements(ydata)
;nfree = n_elements(freearr)-nfixed
;chi2 = chi2 * (npts-n_elements(params) ) / (npts - nfree - 1)

case fail_type of
    0: fail_msg = 'No failure.'
    1: fail_msg = 'Not enough data points.'
    2: fail_msg = 'Loss of precision.'
    3: fail_msg = 'NaN values occurred in calculation.'
    4: fail_msg = 'Too many repeats - unable to step forward.'
    5: fail_msg = 'Maximum number of iterations reached.'
    else: fail_msg = 'Unknown failure code from MCurvefit.'
endcase

; only set # iterations first time we do fit, not after we've refined errors, because
; then it will also be small
if loop_index eq 0 then self ->  set, mcurvefit_iter=iter

self->set, mcurvefit_fail_type=fail_type, $
	mcurvefit_fail_msg=fail_msg, $
	mcurvefit_chi2=chi2, $
	mcurvefit_corr=corr, $
	mcurvefit_covar=covar, $
	mcurvefit_free_mask=freearr

if ~quiet then begin
  print,' '
  ndof = n_elements(ydata) - (n_elements(params) - nfixed) - 1
;  chisq_prob = 100. * chisqr_pdf(chi2*ndof, ndof)
  print,'Mcurvefit results for Interval ' + trim(self.spex_fit_obj->get(/spex_interval_index)) + $
	  ',   Chisq= ' + trim(chi2,'(g12.4)') + ',   Full Chisq= ' + trim(chi2*ndof,'(g12.4)') + $
	  ',   # DoF= ' + trim(ndof) + $ 
;	  ',   # DoF= ' + trim(ndof) + '  Prob.= ' + trim(round(chisq_prob),'(i4)') + '%' + $
    ',   #Iter=' + trim(iter) + $
	  ',   Function ' + self->get(/fit_function,class='fit_function')
  print, 'Params= ', arr2str(trim(params[*],'(g12.4)'), ', ')
  print, 'Sigmas= ', arr2str(trim(sigmas[*],'(g12.4)'), ', ')
  print,' '
endif

end

;------------------------------------------------------------------------------
pro spex_fitalg_mcurvefit__define

struct = { spex_fitalg_mcurvefit, $
           inherits spex_fitalg_gen}
end