;+
; NAME:
;  roboline
; PURPOSE: (one line)
;  Robust fitting of a line
; DESCRIPTION:
;   Robust fitting of y = a = b * x
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  [a,b] = roboline(x, y, sigma, thresh, bplist, GDINDEX = gdindex, coeferr = coefErr)
; INPUTS:
;  x - the indices
;  y - the data array, same length as x
;  sigma - estimated noise for each pixel in y
;  thresh - threshold for outliers 
; OPTIONAL INPUT PARAMETERS:
;  bplist - optional bad pixel list (-1 for all OK)
; INPUT KEYWORD PARAMETERS:
; OUTPUT KEYWORD PARAMETERS:
;  gdindex - named variable that will contain an array of indices used
; OUTPUTS:
;  returns scaleFac such that y ~= a + b * x
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  None.
; PROCEDURE:
;   First iteration finds minimum absolute deviation.
;   Outliers are identified after the first iteration.
;   The second iteration minimizes least squares 
;
; MODIFICATION HISTORY:
;  Written July 22, 2002 Eliot Young, SwRI
;  Rewritten by Leslie A. Young, SwRI, 2002 Oct 29.
;  - first iteration minimizes absolute deviation
;  - exchanged order of x and y in argument list
;  - only two iterations
;  - optionally returns the error and the list of good pixel indices
;-

function roboline, x, y, thresh, sigma=sigma, bplist=bplist, GDINDEX = gdindex, coeferr = coeferr

    n = n_elements(x)
    coefErr = [0.,0.]
    
    if not keyword_set(sigma) then begin
        
        ; First iteration: minimize weighted absolute deviation
        coef = madline(x, y)
        yfit = poly(x,coef)
        
        adev = median(abs(y-yfit))
        sigma = replicate(1.5 * adev, n)
        w = 1/sigma^2
        if keyword_set(bplist) then begin
            if n_elements(bplist) gt 0 then begin
                if bplist[0] ne -1 then w[bplist] = 0.
            end 
        end
        
    endif else begin
  
        ; Get the weights
        w = fltarr(n)
        gd = where(sigma ge 0)
        w[gd] = 1/sigma[gd]^2.
        if keyword_set(bplist) then begin
            if n_elements(bplist) gt 0 then begin
                if bplist[0] ne -1 then w[bplist] = 0.
            end 
        end
        ; First iteration: minimize weighted absolute deviation
        coef = madlinew(x, y, w)
        yfit = poly(x,coef)
        
    endelse
    
    
; Find good points
    gdindex = where( (abs(y-yfit) lt sigma * thresh) and (w ne 0), ngd)
    if ngd eq 0 then begin
        print, 'ROBOLINE: All points rejected in first iteration.'
        return, [0,0]
    end
    
; Second iteration: minimize sum of squares for good points
    sw = total(w[gdindex])
    swx = total(w[gdindex]*x[gdindex])
    swy = total(w[gdindex]*y[gdindex])
    swxy = total(w[gdindex]*x[gdindex]*y[gdindex])
    swxx = total(w[gdindex]*x[gdindex]^2)
    if swxx*sw eq swx^2 then begin
        print, 'ROBOLINEE: indeterminant'
        return, [0,0] 
    end
    det = swxx * sw - swx^2
    a = (swxx*swy - swx * swxy)/det
    b = (-swx*swy + sw * swxy)/det
    aerr = sqrt(swxx/det)
    berr = sqrt(sw/det)
    coeferr = [aerr,berr]
    return, [a,b]
end

pro robolineTEST
    
    n = 101
    x = findgen(n)-(n-1)/2.
    f = x
    a = 0.
    b = 1.
    y = a + b * f
    ysig = replicate(1.,n)
    s = long(systime(1))
    yerr = randomn(s, n) * ysig
    ydat=y+yerr
    for i = 0, 4 do begin
        ibad = floor(randomu(s)*(n-0.0001))
        ydat[ibad] = ydat[ibad]+10*ysig[ibad]*randomu(s)
    end
    thresh = 3.
    bplist = -1
;    robor = roboline(f, ydat, thresh, sigma = ysig, bplist=bplist, GDINDEX = gdindex, coef=roboe)
    robor = roboline(f, ydat, thresh, bplist=bplist, GDINDEX = gdindex, coef=roboe)
    print, a, robor[0],"+/-",roboe[0]
    print, b, robor[1],"+/-",roboe[1]
    plot, x, y, xr=[-.6,.6]*n, yr=[-0.6*n/2.,1.2*max(ydat)], /xs,/ys
    oplot, x, ydat, ps=4
    oplot, x[gdindex], ydat[gdindex], ps=1
    oplot, x, poly(f,robor), line=2, thick=2
end