;+
; NAME:
;  madscalew
; PURPOSE: (one line)
;  minimize the average deviation (MAD) for y = a * x.
; DESCRIPTION:
;  This routine estimates the scaling factor, a, such that
;  the average deviation, total(w * abs(y-a*x) ), is minimized.
;
;  The statistics are robust in that the result is insensitive
;  to outliers.  
;
;  Note that for y[i] = a * x[i], x can be thought of as
;  either a coordinate, or a template (e.g., a line profile).
;
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  a = madscale(x, y)
; INPUTS:
;  x   - dependent variable.
;  y   - data, proportional to x.
; OPTIONAL INPUT PARAMETERS:
;  w   - weights for y
; INPUT KEYWORD PARAMETERS:
; OUTPUT KEYWORD PARAMETERS:
; OUTPUTS:
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  None.
; PROCEDURE:
;
; The goal is to minimize f = total(w * abs(y-a*x)) = total( w * abs(x) * abs(y/x - a) )
; This is equivalent to finding the weighted median of y/x, with
; abs(x) as the weights.
;
; MODIFICATION HISTORY:
;  Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.
;-
function madscalew, x, y, w

  if (n_params() lt 3) then return, madscale(x,y)
  
  gd = where(x ne 0, ngd)
  if ngd eq 0 then return, 0
  return, medianw(y[gd]/x[gd], w * abs(x[gd]))
  
end

pro madscalewTEST

  x = (findgen(21)-10)*0.2
  f = exp(-x^2)
  a = 1000.
  y = a * f
  ysig = sqrt(y)
  yerr = randomn(0L, 21) * ysig
  ydat=y+yerr
  ibad = floor(randomu(0L)*20.99999)
  ydat[ibad] = ydat[ibad]+2*a*randomu(0L)
  mada = madscale(f,ydat)
  madaw = madscalew(f,ydat,1/ysig^2)
  lsa = total(f*ydat)/total(f^2)
  print, a, lsa, mada, madaw
  plot, x, y, xr=[-2.1,2.1], yr=[-0.1,1.2]*max(ydat), /xs,/ys
  oplot, x, ydat, ps=4
  oplot, x, lsa*f, line=1
  oplot, x, madaw*f, line=2, thick=2
end

pro madscalewTEST2
x = [ -2.00, -1.80, -1.60, -1.40, -1.20, -1.00, -0.80, -0.60, -0.40, -0.20,  0.00,  0.20,  0.40,  0.60,  0.80,  1.00,  1.20,  1.40,  1.60,  1.80,  2.00]
f = [  0.02,  0.04,  0.08,  0.14,  0.24,  0.37,  0.53,  0.70,  0.85,  0.96,  1.00,  0.96,  0.85,  0.70,  0.53,  0.37,  0.24,  0.14,  0.08,  0.04,  0.02]
a = 1000.
y = [   18.315639,   39.163883,   77.304726,  140.858429,  236.927750,  367.879456,  527.292358,  697.676331,  852.143738,  960.789429, 1000.000000,  960.789429,  852.143738,  697.676331,  527.292358,  367.879456,  236.927750,  140.858429,   77.304726,   39.163883,   18.315639]
yerr = [    4.741420,    3.161109,  -10.779384,   12.014439,   -0.154639,   24.401287,    0.888818,    1.124579,  -20.697905,   -3.603296,  -29.075243,  -47.393677,   -8.664515,  -33.178108,    4.132473,  -10.847002,  -13.137383,    8.169641,    1.778511,  -10.593134,   -2.154966]
ydat = [   23.057060,   42.324993,   66.525345,  152.872864, 1517.655029,  392.280731,  528.181152,  698.800903,  831.445862,  957.186157,  970.924744,  913.395752,  843.479248,  664.498230,  531.424805,  357.032440,  223.790375,  149.028076,   79.083237,   28.570749,   16.160673]
ysig = sqrt(ydat)


  mada = madscale(f,ydat)
  madaw = madscalew(f,ydat,1/ysig^2)
  lsa = total(f*ydat)/total(f^2)
  print, a, lsa, mada, madaw
  plot, x, y, xr=[-2.1,2.1], yr=[-0.1,1.2]*max(ydat), /xs,/ys
  oplot, x, ydat, ps=4
  oplot, x, lsa*f, line=1
  oplot, x, madaw*f, line=2, thick=2

;xx = ydat/f
;ww = abs(f)/yerr^2
;s = sort(xx)
;xx = xx[s]
;ww = ww[s]/total(ww)

end