;+
; NAME:
;  madline
; PURPOSE: (one line)
;  Fits y = a + bx by the criterion of minimum absolute deviations (MAD).
; DESCRIPTION:
;   Fits y = a + bx by the criterion of least absolute deviations. The arrays x[0..ndata-1] and
;   y[0..ndata-1] are the input experimental points. The fitted parameters a and b are output,
;   along with abdev, which is the mean absolute deviation (in y) of the experimental points from
;   the fitted line. 
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  [a,b] = madline(x, y, yfit, abdev)
; INPUTS:
;  x  - An n-element vector of independent variables.
;  y  - A vector of dependent variables, the same length as X.
; OPTIONAL INPUT PARAMETERS:
; INPUT KEYWORD PARAMETERS:
; OUTPUT KEYWORD PARAMETERS:
; OUTPUTS:
;  returns [a,b], the fitted parameters for y = a + b * x
;  Yfit - A named variable that will contain the vector of calculated Y values. 
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  None.
; PROCEDURE:
;
;
; MODIFICATION HISTORY:
;  Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
;     based on Numerical Recipes medline.
;-

function madline, x, y, yfit, abdev, debug = debug

if not keyword_set(debug) then debug = 0

; find the initial guesses for a and b using least-squares
n = n_elements(x)
sx = total(x)
sy = total(y)
sxx = total(x^2)
sxy = total(x*y)
del = n*sxx-(sx)^2
if del eq 0.0 then begin          ;All X's are the same
    a = median(y, /EVEN) ;Bisect the range w/ a flat line
    b = 0
    yfit = a+b*x
    abdev = mean(abs(y-yfit))
    return, [a,b]
end
a = (sxx*sy-sx*sxy)/del
b = (n*sxy-sx*sy)/del
chisq = total((y-a-b*x)^2)
sigb = sqrt(chisq/del)
if debug then print, 'madline: n, sx, sy, sxx, sxy, del = ',n,sx,sy,sxx,sxy,del
if debug then print, 'madline: initial LS a,b = ',a,b

; get the initial bracket
b1 = b
d=y-b1*x & a = median(d, /EVEN) & f1 = total(x*sgn(d-a))  
b2 = b + 3*sigb*sgn(f1)
d=y-b2*x & a = median(d, /EVEN) & f2 = total(x*sgn(d-a))  
if debug then print, 'madline: initial b = ',b1,b2
  if debug then begin
    yfit = a+b*x
    abdev = mean(abs(y-yfit))
    print, a, b, abdev
  end

; check convergence
if (b1 eq b2) then begin
  yfit = a+b*x
  abdev = mean(abs(y-yfit))
  return, [a,b]
end

; bracket the b's - get f1 and f2 on opposite sides of zero
while ( (f1 * f2) gt 0) do begin
  b = b2 + 1.6 * (b1-b2)
  b1 = b2
  f1 = f2
  b2 = b
  d=y-b2*x & a = median(d, /EVEN) & f2 = total(x*sgn(d-a))  
  if debug then begin
    yfit = a+b*x
    abdev = mean(abs(y-yfit))
    print, a, b, abdev
  end
end
if debug then print, 'madline: bracketing b = ',b1,b2

; refine b via bisection
targ = sigb/100.
while (abs(b2-b1) gt targ) do begin
  if debug then print, 'madline: new bracketing b = ',b1,b2
  b = b1 + 0.5 * (b2-b1)
  if (b eq b1 or b eq b2) then begin
    yfit = a+b*x
    abdev = mean(abs(y-yfit))
    return, [a,b]
  end
  d=y-b*x & a = median(d, /EVEN) & f = total(x*sgn(d-a))  
  if ( (f*f1) ge 0.0) then begin
    f1 = f
    b1 = b
  end else begin
    f2 = f
    b2 = b
  end
  if debug then begin
    yfit = a+b*x
    abdev = mean(abs(y-yfit))
    print, a, b, abdev
  end
end

; converged.  Return.
yfit = a+b*x
abdev = mean(abs(y-yfit))
return, [a,b]

end

pro madlineTEST

  x = [-2.,-1.,0.,1.,2.]
  b = -1
  y = b * x
  yerr = [0.432511,    -0.690602,     -1.03700,      2.23537,     0.505393]*0.2
  ydat=y+yerr
  ydat[1] = 3.5
  res = madline(x,ydat,yfit, abdev)
  print, res
  plot, x, y, xr=[-3,3], yr=[-4,4]
  oplot, x, ydat, ps=4
  oplot, x, yfit, line=2, thick=2
end