; $Id: refinement__define.pro,v 1.5 2006/07/26 14:45:10 azuah Exp $
;###############################################################################
;
; NAME:
;  REFINEMENT__DEFINE
;
; PURPOSE:
;  Object class definition for performing the surface fits in RAINS.
;
; CATEGORY:
;  DAVE, Data Analysis, RAINS, surface fitting
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;+
; NAME:
;       REFINEMENT__DEFINE
;
; PURPOSE:
;
;
; AUTHOR:
;
;       Robert M. Dimeo, Ph.D.
;		NIST Center for Neutron Research
;       100 Bureau Drive
;		Gaithersburg, MD 20899
;       Phone: (301) 975-8135
;       E-mail: robert.dimeo@nist.gov
;       http://www.ncnr.nist.gov/staff/dimeo
;
; CATEGORY:
;
;       Utility, DAVE_REFINE application, DAVE project
;
; CALLING SEQUENCE:
;
;
;
; INPUT PARAMETERS (REQUIRED)
;
;
;
; OPTIONAL OUTPUT KEYWORDS
;
;	ERMSG:		string variable specifying the cause of an error in which the
;				function returns 0.
;
; COMMON BLOCKS:
;
;		None
;
; DEPENDENCIES:
;	MPFIT, MPCURVEFIT, GA_FIT
;
; REQUIREMENTS:
;
;	Written using IDL 5.6.  I cannot ensure that it will work with previous
;	releases.
;
; METHODS:
;
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, April 9, 2003
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro refinement::cleanup
ptr_free,self.pfit_info,self.pga_info,self.res_ptr
ptr_free,self.iterargs_ptr,self.pcor_ptr
ptr_free,self.yindex_ptr,self.xindex_ptr,self.zfit_ptr,self.residuals_ptr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro refine_fitfunction,x,parms,yfit,_Extra = extra
omodel = extra.omodel & odata = extra.odata
yindex = *extra.yindex_ptr
xindex = *extra.xindex_ptr
omodel->setparms,parms
ret = omodel->evaluate(zvalues = zvalues)
yfit = (zvalues[xindex,*])[*,yindex]
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;function refinement::set_ga_parms,	funcname,							$
;									PARAM_BOUNDS = param_bounds,		$
;									NGEN = ngen,						$
;									NITER = niter,						$
;									KEEP_BEST = keep_best,				$
;									NPOP = npop,						$
;									PCROSS = pcross,					$
;									PMUT = pmut,						$
;									AMP_MUT = amp_mut,					$
;									ERRMSG = ermsg,						$
;									PARMS = parms,						$
;									AVE_FIT = ave_fit,					$
;									BEST_FIT = best_fit,				$
;									GA_BEST_FIT = ga_best_fit,			$
;									LS_FIT = ls_fit,					$
;									SIGMA = sigma
;
;return,1
;end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::set_ls_parms,	xtol = xtol,	$
									gtol = gtol,	$
									ftol = ftol,	$
									itmax = itmax
if n_elements(xtol) ne 0 then self.xtol = xtol
if n_elements(gtol) ne 0 then self.gtol = gtol
if n_elements(ftol) ne 0 then self.ftol = ftol
if n_elements(itmax) ne 0 then self.itmax = itmax
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::get_ls_parms,	xtol = xtol,	$
									gtol = gtol,	$
									ftol = ftol,	$
									itmax = itmax
xtol = self.xtol
gtol = self.gtol
ftol = self.ftol
itmax = self.itmax
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::set_ga_parms
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::set_fit_parms,errmsg = errmsg
ret = 1 & errmsg = 'No errors generated in setting the fit parameters'
case self.algorithm of
'LS':	$
	begin
		ret = self->set_ls_parms()
		if ret eq 0 then begin
			errmsg = 'Least-squares parameters could not be set'
		endif
	end
'GA':	$
	begin
		ret = self->set_ga_parms()
		if ret eq 0 then begin
			errmsg = 'Genetic algorithm parameters could not be set'
		endif
	end
'GA_LS':	$
	begin
		ret = self->set_ga_parms()
		if ret eq 0 then begin
			errmsg = 'Genetic algorithm parameters could not be set'
			return,0
		endif
		ret = self->set_ls_parms()
		if ret eq 0 then begin
			errmsg = 'Least-squares parameters could not be set'
			return,0
		endif
	end
else:
endcase
return,ret
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::display_residuals
; Plot the residuals
if n_elements(*self.residuals_ptr) eq 0 then return,0
ret = self.odata->get(nx = nx,ny = ny,xvals = x,yvals = y,qty = qty,err = err)
xrange = [min(x),max(x)]
yrange = [min(y),max(y)]
if n_elements(*self.residuals_ptr) eq 0 then return,0
contour,bytscl(*self.residuals_ptr),x[*self.xindex_ptr],y[*self.yindex_ptr], $
	/fill,nlevels = 20, $
	xrange = xrange, /xsty, $
	yrange = yrange,/ysty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::display_image
ret = self.odata->get(nx = nx,ny = ny,xvals = x,yvals = y,qty = qty,err = err)
xrange = [min(x),max(x)]
yrange = [min(y),max(y)]
; display the data then overlay the model
contour,bytscl(alog10(0.1+qty)),x,y,/fill,nlevels = 20, $
	xrange = xrange, /xsty, $
	yrange = yrange,/ysty
if n_elements(*self.zfit_ptr) ne 0 then $
	contour,bytscl(alog10(0.1+*self.zfit_ptr)), $
	x[*self.xindex_ptr],y[*self.yindex_ptr],nlevels = 15,/overplot

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function status_message,status
if n_params() eq 0 then return,0
case status of
0:  status_string = 'improper input parameters'
1:	status_string = 'both actual and predicted relative reductions' $
					+' in the sum of squares are at most FTOL'
2:	status_string = 'relative error between two consecutive iterates' $
					+' is at most XTOL'
3:	status_string = 'conditions for STATUS = 1 and STATUS = 2 both hold'
4:	status_string = 'the cosine of the angle between fvec and any' $
					+' column of the jacobian is at most GTOL in absolute value'
5:	status_string = 'the maximum number of iterations has been reached'
6:	status_string = 'FTOL is too small. no further reduction in' $
					+' the sum of squares is possible'
7:	status_string = ' XTOL is too small. no further improvement in' $
					+' the approximate solution x is possible'
8:	status_string = ' GTOL is too small. fvec is orthogonal to the' $
					+' columns of the jacobian to machine precision'
else:	status_string = ''
endcase
return,status_string
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::fit_ga
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro refine_iterproc,fnc,p,iter,fnorm,				$
			stopBut = stopBut,txt_field =txt_field, $
			chisq_field = chisq_field,	$
			FUNCTARGS = functargs,		$
			PARINFO =  parinfo,			$
			QUIET = quiet,				$
			_Extra = extra

if n_elements(stopBut) eq 0 then return
if widget_info(stopBut,/valid_id) eq 0 then return
widget_control,txt_field,set_value = strtrim(string(iter),2)
widget_control,chisq_field,set_value = strtrim(string(fnorm/extra.dof),2)
event = widget_event(/nowait,stopBut)
evname = tag_names(event,/structure_name)
if evname eq '' or event.id eq 0 then return

if evname eq 'WIDGET_BUTTON' and event.select then begin
  event1 = widget_event(/nowait,stopBut)
  common mpfit_error,mperr
  mperr = -1
endif

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::fit_ls
ret = self.odata->get(nx = nx,ny = ny,xvals = x,yvals = y,qty = qty,err = err)
data=(qty[*self.xindex_ptr,*])[*,*self.yindex_ptr]
error=(err[*self.xindex_ptr,*])[*,*self.yindex_ptr]
x = x[*self.xindex_ptr]
y = y[*self.yindex_ptr]
weights = 1.d/error^2
parms = self.omodel->getparms()
parinfo = self.omodel->createparinfo()
tstart = systime(/seconds)
yf = mpcurvefit([x,y], double(data),double(weights), parms, sigma, $
                      FUNCTION_NAME='refine_fitfunction', $
                      status = status,								$
                      /quiet, 										$
                      chisq = chisq, 								$
                      parinfo = parinfo, 							$
                      itmax = self.itmax,							$
                      ftol = self.ftol,								$
                      gtol = self.gtol,								$
                      xtol = self.xtol,								$
                      covar = covar,								$
                      iterargs = *self.iterargs_ptr,				$
                      iterproc = 'refine_iterproc',					$
                      functargs = {	odata:self.odata,				$
                      				omodel:self.omodel,				$
                      				xindex_ptr:self.xindex_ptr,		$
                      				yindex_ptr:self.yindex_ptr	},	$
                      /autoderivative, $
                      ERRMSG=errmsg)
if errmsg ne '' then return,0
; Compute the correlation matrix
PCOR = COVAR * 0
FOR i = 0, n_elements(parms)-1 DO FOR j = 0, n_elements(parms)-1 DO $
	PCOR(i,j) = COVAR(i,j)/sqrt(COVAR(i,i)*COVAR(j,j))
*self.pcor_ptr = pcor

dof = n_elements(data) - n_elements(parms)
self.chisq = chisq/dof
self.omodel->setchisq,self.chisq
self.omodel->setparmError,sigma
*self.zfit_ptr = yf
*self.residuals_ptr = (data-yf)/err
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::fit_report,output = output
self.omodel->displayparms,output = output
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::fit_gals
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::fit
case self.algorithm of
	'LS':		ret = self->fit_ls()
	'GA':		ret = self->fit_ga()
	'GA/LS':	ret = self->fit_gals()
	else:
endcase
return,ret
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::set,	resolution = resolution,		$
							algorithm = algorithm,			$
							xindex = xindex,				$
							yindex = yindex,				$
							iterproc = iterproc,			$
							iterargs = iterargs
if n_elements(iterproc) ne 0 then self.iterproc = iterproc
if n_elements(iterargs) ne 0 then *self.iterargs_ptr = iterargs
if n_elements(xindex) ne 0 then *self.xindex_ptr = xindex
if n_elements(yindex) ne 0 then *self.yindex_ptr = yindex
if n_elements(resolution) ne 0 then *self.res_ptr = resolution
if n_elements(algorithm) ne 0 then self.algorithm = algorithm
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::get,	algorithm = algorithm,		$
							chisq = chisq,				$
							pcor = pcor,				$
							xindex = xindex,			$
							yindex = yindex

xindex = *self.xindex_ptr
yindex = *self.yindex_ptr
if arg_present(algorithm) then algorithm = self.algorithm
if arg_present(chisq) then chisq = self.chisq
if arg_present(pcor) then begin
	if n_elements(*self.pcor_ptr) gt 0 then begin
		pcor = *self.pcor_ptr
	endif else begin
		return,0
	endelse
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function refinement::init,	omodel = omodel,			$
							odata = odata,				$
							resolution = resolution,	$
							xindex = xindex,			$
							yindex = yindex,			$
							iterproc = iterproc,		$
							iterargs = iterargs,		$
							algorithm = algorithm,		$
							xtol = xtol,				$
							itmax = itmax,				$
							gtol = gtol,				$
							ftol = ftol,				$
							errmsg = errmsg
errmsg = 'No errors generated in refinement initialization'
if n_elements(algorithm) eq 0 then algorithm = 'LS'
if n_elements(omodel) eq 0 then begin
	errmsg = 'You must pass in a model object'
	return,0
endif
if n_elements(odata) eq 0 then begin
	errmsg = 'You must pass in a data object'
	return,0
endif
if (not obj_valid(omodel)) then begin
	errmsg = 'You must pass in a valid model object'
	return,0
endif
if (not obj_valid(odata)) then begin
	errmsg = 'You must pass in a valid data'
        return,0
     endif
ret = odata->get(nx = nx,ny = ny,xvals = x,yvals = y,qty = qty,err = err)

self.iterargs_ptr = ptr_new(/allocate_heap)
if n_elements(iterargs) ne 0 then *self.iterargs_ptr = iterargs
if n_elements(iterproc) ne 0 then self.iterproc = iterproc
if n_elements(xindex) eq 0 then xindex = indgen(nx)
if n_elements(yindex) eq 0 then yindex = indgen(ny)
if n_elements(xindex) lt 2 or n_elements(yindex) lt 2 then begin
   errmsg = "There aren't enough points to fit in the selected range"
   return,0
endif
self.odata = odata
self.omodel = omodel
self.xindex_ptr = ptr_new(xindex,/no_copy)
self.yindex_ptr = ptr_new(yindex,/no_copy)
self.pfit_info = ptr_new(/allocate_heap)
self.pcor_ptr = ptr_new(/allocate_heap)

self.pga_info = ptr_new(/allocate_heap)
self.res_ptr = ptr_new(/allocate_heap)
self.zfit_ptr = ptr_new(/allocate_heap)
self.residuals_ptr = ptr_new(/allocate_heap)
self.chisq = 0.0
if n_elements(itmax) eq 0 then itmax = 200
if n_elements(ftol) eq 0 then ftol = 1D-10
if n_elements(gtol) eq 0 then gtol = 1D-10
if n_elements(xtol) eq 0 then xtol = 1D-10
self.itmax = itmax
self.ftol = ftol
self.gtol = gtol
self.xtol = xtol
if n_elements(resolution) ne 0 then begin
   *self.res_ptr = resolution
endif
self.algorithm = algorithm
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro refinement__define
  define =	{	refinement,					$
                        algorithm:'',				$
                        chisq:0.0,					$
                        itmax:0,					$
                        ftol:0D,					$
                        gtol:0D,					$
                        xtol:0D,					$
                        pcor_ptr:ptr_new(),			$
                        pfit_info:ptr_new(),		$
                        pga_info:ptr_new(),			$ ; holds the GA parameter info
                        omodel:obj_new(),			$
                        res_ptr:ptr_new(),			$
                        zfit_ptr:ptr_new(),			$
                        residuals_ptr:ptr_new(),	$
                        yindex_ptr:ptr_new(),		$
                        xindex_ptr:ptr_new(),		$
                        iterproc:'',				$
                        iterargs_ptr:ptr_new(),		$
                        odata:obj_new()				}
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
