;+
; NAME:
;  madlinew
; 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,
;   which minimize the sum of (w * |y - a - bx| ). 
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  [a,b] = madlinew(x, y, w, yfit)
; INPUTS:
;  x  - An n-element vector of independent variables.
;  y  - A vector of dependent variables, the same length as X.
;  w  - weights
; 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 medfit.
;-

function madlinew, x, y, w, yfit, abdev

debug = 0

; find the initial guesses for a and b using least-squares
sw = total(w)
swx = total(w*x)
swy = total(w*y)
swxx = total(w*x^2)
swxy = total(w*x*y)
del = sw*swxx-(swx)^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
    return, [a,b]
end
a = (swxx*swy-swx*swxy)/del
b = (sw*swxy-swx*swy)/del
chisq = total((y-a-b*x)^2*w)
sigb = sqrt(chisq/del)
if debug then print, 'MADLINEW: initial LS a,b = ',a,b

; get the initial bracket
b1 = b
d=y-b1*x & a = medianw(d,w) & f1 = total(w*x*sgn(d-a))  
b2 = b + 3*sigb*sgn(f1)
d=y-b2*x & a = medianw(d,w) & f2 = total(w*x*sgn(d-a))  
if debug then print, 'MADLINEW: initial b = ',b1,b2

; 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 = medianw(d,w) & f2 = total(w*x*sgn(d-a))  
end
if debug then print, 'MADLINEW: bracketing b = ',b1,b2

; refine b via bisection
targ = sigb/100.
while (abs(b2-b1) gt targ) do begin
  if debug then print, 'MADLINEW: 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 = medianw(d,w) & f = total(w*x*sgn(d-a))  
  if ( (f*f1) ge 0.0) then begin
    f1 = f
    b1 = b
  end else begin
    f2 = f
    b2 = b
  end
end

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

end

pro madlinewTEST

  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
  w = [1,1,1,1,1]*2
  res = madlinew(x,ydat,w, 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


pro madlinewTEST2

  x1 = findgen(10)
  x2 = findgen(10) + 20
  x3 = findgen(10) + 40
  x = [x1,x2,x3]
  b = -1
  y = b * x
  ysig = [replicate(0.1,10), replicate(0.1,10), replicate(2,10)]
  yerr = randomn(0L, 30)*ysig
  ydat=y+yerr
  ydat[1] = 3.5
  w = 1/ysig^2
  res = madlinew(x,ydat,w, yfit, abdev)
  print, res
  plot, x, y, xr=[-2, 52], yr=[-60,10], /xs, /ys
  oplot, x, ydat, ps=4
  oplot, x, yfit, line=2, thick=2
end