pro Max_Resid_Like, image_deconv, chi_sq, conv, LAG_AUTO_CORR=Lag_ac,  $
						SIGMA_NOISE=sig_noise, $
						INITIALIZE=init,       $
						IMAGE_DATA=image_data, $
						PSF_DATA=psf_data,     $
						POSITIVITY_EPS=eps
;+
; NAME:
;	Max_Resid_Like
; PURPOSE:
; CALLING EXAMPLE:
; INPUTS:
; KEYWORDS:
; OUTPUTS:
; EXTERNAL CALLS:
; COMMON BLOCKS:
; PROCEDURE:
;	Uses conjugate-gradient method to minimize
;	the chi-square of residual auto-correlation,
;	as computed by function chi_Resid_Corr( image_deconv, gradient ).
;	The chi_sq is minimized until close to most probable value (mode),
;	and then the difference (chi_sq - mode)^2 is further minimized
;	(phase2 of algorithm).
; Based on:
;	"Incorporation of Spatial Information in Bayesian Image Reconstruction:
;			the Maximum Residual Likelihood (MRL) Criterion",
;	by Robert Pina and Richard Puetter UCSD 1992.
; MODIFICATION HISTORY:
;	Written, Frank Varosi NASA/GSFC 1992.
;-
  common Deconv_Data, imaged, psf, psf_FT
  common chi_Resid_Corr, sigma, Lagac, mode
  common Max_Resid_Like, close, phase2
  common chi_Resid_Corr1, epsilon

	if N_elements( eps ) EQ 1 then epsilon=eps else epsilon=1.e-7
	if N_elements( image_data ) GT 1 then imaged = image_data

	if N_elements( psf_data ) GT 1 then begin
		psf = psf_data
		psf_FT = 0
	   endif

	if keyword_set( init ) then begin
		close = 0
		phase2 = 0
		if N_elements( sig_noise ) GE 1 then sigma = sig_noise $
						else sigma = 0
		if N_elements( Lag_ac ) GT 0 then Lagac = Lag_ac $
			else Lagac = round_off( FullWid_HalfMax( psf,/AV ) ) > 1
		sim = size( imaged )
		Lagac = Lagac < ( min( sim(1:2) )/2 -1 )
	   endif

	if N_elements( image_deconv ) NE N_elements( imaged ) then begin
		init = 1
		image_deconv = imaged
		image_deconv(*) = total( imaged )/N_elements( imaged )
	   endif

	if (close) AND (NOT phase2) then begin
		phase2 = 1
		init = 1
		message,"MRL is close to convergence, entering phase II",/INFO
	   endif

	minF_conj_grad, image_deconv, chi_sq, conv, FUNC="chi_Resid_Corr", $
						INIT=init, TOLER=1.e-3,    $
						USE_DERIV=(epsilon GT 0)

	if (epsilon LE 0) then image_deconv = image_deconv > 0
	close = chi_sq  LT  ( mode + sqrt( 2*(mode+2) ) )
return
end