; $Id$
;
;+
; NAME:
;       CG_OPTIMIZE
;
; LAST CHANGED: ------------------------------  28-jul-01 (aph)
;
; PURPOSE:
;       This function finds the minimum of a user-specified objective
;       function using the conjugate gradient algorithm (Polak-Ribiere
;       variant). The objective function to be minimized must return a
;       scalar-valued output given a vector-valued input. Also required
;       is a function to calculate the gradient of the objective function.
;
; CATEGORY:
;       Optimization
;
; CALLING SEQUENCE:
;       x = CG_OPTIMIZE(xstart, Func, DFunc)
;
; INPUTS:
;       xstart:    A starting guess for the solution vector x. A common
;                  starting guess is the zero vector.
;
;       Func:      A string containing the name of the objective function to
;                  minimize. The objective function must take a vector-valued
;                  input and return a scalar-valued output.
;
;       DFunc:     A string containing the name of a function that calculates
;                  the gradient of Func at a vector-valued input point. The
;                  output of the function must be a vector of the same length
;                  as the input vector.
;
; KEYWORD PARAMETERS:
;       NUM_ITER:  The maximum number of iterations of the algorithm to run.
;                  The default is 100.
;
;       TOLERANCE: If set, iterates will be computed until the change in
;                  the objective function is smaller than the tolerance,
;                  or the maximum number of iterations is reached, whichever
;                  comes first.
;
;       OBJECTIVE: A variable to receive the value of the objective function
;                  at each iteration. It will therefore be a vector of length
;                  one plus the number of iterations performed.
;
;       LINMIN:    A string containing the name of a function that performs
;                  one dimensional minimization of an objective function along
;                  a vector direction from a starting point. The function
;                  must take the following arguments:
;                    (start, dir, fmin, Func, DFunc, _EXTRA)
;                  where start is the starting point (vector), dir is the
;                  direction (vector), fmin is the value of the objective
;                  function at the line minimum (scalar), Func, DFunc,
;                  and _EXTRA are as described for CG_OPTIMIZE.
;                  The 1-D minimization function must return a vector that
;                  minimizes the objective function along the specified
;                  direction.
;                  If not specified, the default function CG_LINMIN is used.
;
;		POS_ONLY:  If set, force all values to be optimized to be GT 0
;
;		LIMITS:		[x-min, x-max] lower and upper bounds
;
;		SWAP:		If set, replace low values with upper limit and vice-verso
;
;       _EXTRA:    Keywords to Func and DFunc may be passed to those functions
;                  through the CG_OPTIMIZE function. This is useful for
;                  providing data to the objective and gradient functions.
;
; OUTPUTS:
;       x:         The solution vector that minimizes Func.
;
; ROUTINES INCLUDED IN THIS MODULE:
;       CG_LINMIN   The default line minimization routine.
;       CG_OPTIMIZE The main function.
;
;       and example functions: CGEX_FUNC, CGEX_DFUNC, CGEX_LINMIN
;
; WRITTEN BY:
;       Billy W. Loo, Jr.
;       Bioengineering Graduate Group, UCSF / UCB
;       School of Medicine, UCD
;       Lawrence Berkeley National Lab
;       February, 2000
;
; MODIFICATION HISTORY:
; (29-dec-00 aph) add ALLPOS to force non-negative values (not working as parameter)
; (13-feb-01 aph) adapt for use in a curve-fit routine
; (27-jul-01 aph) get pos_only to work as an option; - executes but no convergence
;              explicitly include A, B in all call functions
; (28-jul-01 aph) add upper / lower bound constraints and swap
;-

; ***************************************

FUNCTION CGEX_LINMIN, x, v, fmin, Func, DFunc, A=a, B=b
  d  = A ## x - b
  Av = A ## v

  lambda = TOTAL(d * Av) / TOTAL(Av * Av)

  ;Update difference vector.
  d = TEMPORARY(d) - lambda * Av

  fmin = TOTAL(d * d)

  RETURN, x - lambda * v
END

; ***************************************

;Define objective function using least squares measure.
FUNCTION CGEX_FUNC, x, A=a, B=b
  d = A ## x - b

  RETURN, TOTAL(d * d)
END

; ***************************************

FUNCTION CGEX_DFUNC, x, A=a, B=b
  d = A ## x - b

  RETURN, TRANSPOSE(A) ## d
END

; ***************************************

FUNCTION CG_OPTIMIZE, xstart, Func, DFunc, $
                      NUM_ITER=num_iter,   $
                      TOLERANCE=tolerance, $
                      START=start,         $
                      OBJECTIVE=objective, $
                      LINMIN=Linmin,       $
                      A = a,               $
                      B = b,               $
                      POS_only=pos_only,   $
                      LIMITS = limits,     $
                      SWAP = swap,         $
                      _EXTRA=extra

  ;Set default number of iterations.
  IF (N_ELEMENTS(num_iter) EQ 0) THEN num_iter = 100

  ;Set default tolerance.
  IF (N_ELEMENTS(tolerance) EQ 0) THEN BEGIN
    tolerance = MACHAR()
    tolerance = SQRT(tolerance.eps)
  ENDIF

  ;Set default line minimization function.
  IF (N_ELEMENTS(linmin) EQ 0) THEN Linmin = 'CG_LINMIN'

  ;Set starting solution vector.
  x = xstart

  ;Calculate the starting objective function value.
  objective    = FLTARR(num_iter + 1)
  objective[0] = CALL_FUNCTION(Func, x, A=a,b=b, _EXTRA=extra)

  ;Calculate starting gradient.
  g = CALL_FUNCTION(DFunc, x, A=a, b=b, _EXTRA=extra)

  ;Starting search direction.
  v = -g

  ;Iterate over conjugate search directions.
  FOR i = 0L, num_iter-1 DO BEGIN

; ------ force x to be positive or zero, if requested ---- THIS DIVERGES !!
	if keyword_set(pos_only) then begin
		FOR ii = 0, n_elements(x)-1 do begin
	   		if x(ii) LT 0 then begin
;	   		   print, ii, x(ii), ' set to zero'
	   		   	x(ii) = 0.
	   		endif
	   ENDFOR
	endif

; ------ enforce limits if requested ------
	if keyword_set(limits) then begin
		if NOT keyword_set(swap) then begin
			x = min([x,limits(0)])
			x = max([x,limits(1)])
		endif else begin
			FOR ii = 0, n_elements(x)-1 do begin
		   		if x(ii) LT limits(0) then x(ii) = limits(0)
		   		if x(ii) GT limits(1) then x(ii) = limits(1)
	   		ENDFOR
	   	endelse
	endif

; Minimize along search direction.
;    print, 'first time in cg_opt'  & help, x, v
    x = CALL_FUNCTION(Linmin, x, v, fmin, Func, DFunc, A=a, b=b, _EXTRA=extra)
;    help, x,v

;Update the objective function value.
    objective[i+1] = fmin

;Quit if tolerance has been reached.
    IF (ABS(objective[i+1] - objective[i]) LE tolerance) THEN BEGIN
	      objective = objective[0:i+1]
	      GOTO, Done
    ENDIF

;New gradient.
    ng = CALL_FUNCTION(DFunc, x,  A=a, b=b,_EXTRA=extra)

;Scaling factor for search direction computation.
    ;(Polak-Ribiere version of CG algorithm)
    Gamma = TOTAL((ng - g) * ng) / TOTAL(g * g)

;Update gradient.
    g = TEMPORARY(ng)

;Update search direction.
    v = Gamma * v - g
  ENDFOR

  Done:

  num_iter = N_ELEMENTS(objective) - 1

  RETURN, x
END ;cg_optimize
