Function amofitcall, func, p, coords, data, weights, method, model
; evaluates the function and the quantity to be minimized

  call_procedure,func, coords, p, model  ;Eval fcn at trial point
  r = data-model; residual
  case method of 
    'maxabres': y = max(abs(r))
    'sumabres': y = total(weights*abs(r))
    'sumsqres': y = total(weights*r^2)
    else: y = total(weights*r^2)  ; default to the old standby, sum of squared residuals
  endcase
  return, y
end


Function amofittry, p, y, psum, func, ihi, fac, coords, data, weights, method, model
; Extrapolates by a factor fac through the face of the simplex, across
; from the high point, tries it and replaces the high point if the new
; point is better.

  fac1 = (1.0 - fac) / n_elements(psum)
  fac2 = fac1  - fac
  ptry = psum * fac1 - p[*,ihi] * fac2
;  ytry = call_function(func, ptry)  ;Eval fcn at trial point
  ytry = amofitcall(func, ptry, coords, data, weights, method, model)  ;Eval fcn at trial point
  if ytry lt y[ihi] then begin    ;If its better than highest, replace highest
    y[ihi] = ytry
    psum = psum + ptry - p[*,ihi]
    p[0,ihi] = ptry
    endif
  return, ytry
end


Function Amoebafit, ftol, data, COORDS=coords, WEIGHTS=weights, METHOD=method, FUNCTION_NAME=func, FUNCTION_VALUE=y, $
	NCALLS = ncalls, NMAX = nmax, P0 = p0, SCALE=scale, SIMPLEX=p, MODEL_VALUE = model
; The Numerical Recipes procedure Amoeba, with some embellishments.
;
; Reference: Numerical Recipes, 2nd Edition, Page 411.
;  P = fltarr(ndim, ndim+1)
;+
; NAME:
;	AMOEBA
;
; PURPOSE:
;	Multidimensional minimization of a function FUNC(X), where
;	X is an N-dimensional vector, using the downhill simplex
;	method of Nelder and Mead, 1965, Computer Journal, Vol 7, pp 308-313.
;
;    This routine is modified from IDL's amoeba routine, which
;	in turn is based on the AMOEBA routine, Numerical
;	Recipes in C: The Art of Scientific Computing (Second Edition), Page
;	411.
;
; CATEGORY:
;	Function minimization/maximization. Simplex method.
;
; CALLING SEQUENCE:
;	Result = AMOEBAFIT(Ftol, data, ....)
; INPUTS:
;    FTOL:  the fractional tolerance to be achieved in the function
;	value.  e.g. the fractional decrease in the function value in the
;	terminating step.  This should never be less than the
;	machine's single or double precision.
;    DATA: The data values
; KEYWORD PARAMETERS:
;    FUNCTION_NAME: a string containing the name of the function to
;	be minimized.  If omitted, the function FUNC is minimized.
;	This function must accept an Ndim vector as its only parameter and
;	return a scalar single or double precision floating point value as its
;	result. 
;    *************************************************************
;    This function is now the same form as for curvefit, eg
;    function_name(coords, param, values, deriv)
;    *************************************************************
;
;    COORD: The dependent variable, [0,1,2,3...] if not passed
;    WEIGHT: The weights for the data, [1,1,1,1,] if not passed
;
;    FUNCTION_VALUE: (output) on exit, an Ndim+1 element vector
;	containing the function values at the simplex points.  The first
;	element contains the function minimum. 
;    
;    NCALLS: (output) the of times the function was evaluated. 
;    NMAX: the maximum number of function evaluations allowed
;	before terminating.  Default = 5000.
;    P0: Initial starting point, an Ndim element vector.  The starting
;	point must be specified using either the keyword SIMPLEX, or P0 and
;	SCALE.  P0 may be either single or double precision floating.
;	For example, in a 3-dimensional problem, if the initial guess
;	is the point [0,0,0], and it is known that the function's
;	minimum value occurs in the interval: -10 <
;	X(0) < 10, -100 < X(1) < 100, -200 < X(2) < 200, specify: P0=[0,0,0],
;	SCALE=[10, 100, 200]. 
;    SCALE: a scalar or Ndim element vector contaiing the problem's
;	characteristic length scale for each dimension.
;	SCALE is used with P0 to form an initial (Ndim+1) point simplex.
;	If all dimensions have the same	scale, specify a scalar.
;    SIMPLEX: (output and/or optional input) On input, if P0 and SCALE
;	are not set, SIMPLEX contains the Ndim+1 vertices, each of
;	Ndim elements, of starting simplex, in either single or double
;	precision floating point, in an (Ndim, Ndim+1) array. On output,
;	SIMPLEX contains the simplex, of dimensions (Ndim, Ndim+1), enclosing
;	the function minimum.  The first point, Simplex(*,0), corresponds to
;	the function's minimum.
;	METHOD: What to minimize.  Choices are 
;       'sumsqres' to minimize the sum of (data-model)^2, 
;       'sumabres' to minimize the sum of |data-model|, 
;       'maxabdev' to minimize the maximum of |data-model|
;
; OUTPUTS:
;   Result: If the minimum is found, an Ndim vector, corresponding to
;	the Function's minimum value is returned.  If a function minimum
;	within the given tolerance, is NOT found in the given number of
;	evaluations, a scalar value of -1 is returned.
;
; COMMON BLOCKS:
;	None.
;
; SIDE EFFECTS:
;	None.
;
; PROCEDURE:
;	This procedure implements the Simplex method, described in
;	Numerical Recipes, Section 10.4.  See also the POWELL procedure.
;
;	Advantages:  requires only function evaluations, not
;	derivatives, may be more reliable than the POWELL method.
;	Disadvantages: not as efficient as Powell's method, and usually
;	requires more function evaluations.
;
;	Results are performed in the mode (single or double precision)
;	returned by the user-supplied function.  The mode of the inputs P0,
;	SCALE, or SIMPLEX, should match that returned by the function. The
;	mode of the input vector supplied to the user-written function, is
;	determined by P0, SCALE, or SIMPLEX.
;
; EXAMPLE:
;	Use Amoeba to find the slope and intercept of a straight line fitting
;	a given set of points minimizing the maximum error:
;
;	The function to be minimized returns the maximum error,
;	given p(0) = intercept, and p(1) = slope:
; FUNCTION FUNC, x, p, y
; RETURN, (p[0] + p[1] * x)
; END
;
;	Define the data points:
;   x = findgen(17)*5
;   y = [ 12.0,  24.3,  39.6,  51.0,  66.5,  78.4,  92.7, 107.8, 120.0, $
;        135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]
;
;	Call the function.  Fractional tolerance = 1 part in 10^5, 
;	Initial guess = [0,0], and the minimum should be found within
;	a distance of 100 of that point: 
;   r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_VALUE=fval, METHOD='maxabres')
;
;	Check for convergence:
;   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
;	Print results.
;   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
;Intercept, Slope:      11.4100      2.72800
;Function value:       1.33000
;

;
; MODIFICATION HISTORY:
;	LAY, Oct, 2002.	Modified Amoeba.
;-

if keyword_set(scale) then begin  ;If set, then p0 is initial starting pnt
   ndim = n_elements(p0)
   p = p0 # replicate(1.0, ndim+1)
   for i=0, ndim-1 do p[i,i+1] = p0[i] + scale[i < (n_elements(scale)-1)]
endif

s = size(p)
if s[0] ne 2 then message, 'Either (SCALE,P0) or SIMPLEX must be initialized'
ndim = s[1]			;Dimensionality of simplex
mpts = ndim+1			;# of points in simplex
;if n_elements(func) eq 0 then func = 'FUNC'
if n_elements(func) eq 0 then func = 'FUNCT'  ; LAY default function for curvefit
if n_elements(nmax) eq 0 then nmax = 5000L
; Start LAY changes
ndata = n_elements(data) 
if n_elements(coords) ne ndata then coords = findgen(ndata)
if n_elements(weights) ne ndata then weights = replicate(1., ndata)
if n_elements(method) eq 0 then method = 'sumsqres'

; End LAY changes

;y = replicate(call_function(func, p[*,0]), mpts)  ;Init Y to proper type
;for i=1, ndim do y[i] = call_function(func, p[*,i]) ;Fill in rest of the vals
; Start LAY changes
;Init Y to proper type. LAY new function call
y = replicate(amofitcall(func, p[*,0], coords, data, weights, method), mpts)
for i=1, ndim do y[i] = amofitcall(func, p[*,i], coords, data, weights, method) ;Fill in rest of the vals
; End LAY changes
ncalls = 0L
psum = total(p,2)

while ncalls le nmax do begin		;Each iteration

  s = sort(y)
  ilo = s[0]		;Lowest point
  ihi = s[ndim]		;Highest point
  inhi = s[ndim-1]	;Next highest point
  d = abs(y[ihi]) + abs(y[ilo]) ;Denominator = interval
  if d ne 0.0 then rtol = 2.0 * abs(y[ihi]-y[ilo])/d $
  else rtol = ftol / 2.         ;Terminate if interval is 0

  if rtol lt ftol then begin ;Done?
      t = y[0] & y[0] = y[ilo] & y[ilo] = t ;Sort so fcn min is 0th elem
      t = p[*,ilo] & p[*,ilo] = p[*,0] & p[*,0] = t
      temp = amofitcall(func, p[*,ilo], coords, data, weights, method, model) ; fill in model
      return, t                 ;params for fcn min
  endif
    
  ncalls = ncalls + 2
  ytry = amofittry(p, y, psum, func, ihi, -1.0, coords, data, weights, method, modeltry)
  if ytry le y[ilo] then ytry = amofittry(p,y,psum, func, ihi, 2.0, coords, data, weights, method, modeltry) $
  else if ytry ge y[inhi] then begin
    ysave = y[ihi]
    ytry = amofittry(p,y,psum,func, ihi, 0.5, coords, data, weights, method, modeltry)
    if ytry ge ysave then begin
	for i=0, ndim do if i ne ilo then begin
	  psum = 0.5 * (p[*,i] + p[*,ilo])
	  p[*,i] = psum
;	  y[i] = call_function(func, psum)
	  y[i] = amofitcall(func, psum, coords, data, weights, method)
	  endif
	ncalls = ncalls + ndim
	psum = total(p,2)
	endif		;ytry ge ysave
    endif else ncalls = ncalls  - 1
  endwhile
  return, -1		;Here, the function failed to converge.
end

pro amofitfunc, x, p, y
 y = (p[0] + p[1] * x)
END

pro amoebafittest
;	Define the data points:
   x = findgen(17)*5
   y = [ 12.0,  24.3,  39.6,  51.0,  66.5,  78.4,  92.7, 107.8, 120.0, $
        135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]
;
;	Call the function.  Fractional tolerance = 1 part in 10^5, 
;	Initial guess = [0,0], and the minimum should be found within
;	a distance of 100 of that point: 
   r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_NAME='amofitfunc', FUNCTION_VALUE=fval, METHOD='maxabres')
;
;	Check for convergence:
   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
;	Print results.
   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
;Intercept, Slope:      11.4100      2.72800
;Function value:       1.33000


   r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_NAME='amofitfunc', FUNCTION_VALUE=fval, METHOD='sumabres')
   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)

   r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_NAME='amofitfunc', FUNCTION_VALUE=fval, METHOD='sumsqres')
   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
end