;+
; PROJECT:
;   HESSI
;
; NAME:
;   hsi_conditional_average
;
; PURPOSE:
;   Returns the average value of an nv-element input array, excluding values that differ from final average by more than limitfactor*average
;
; CATEGORY:
;   HESSI/IDL/UTIL;
;
; CALLING SEQUENCE:
;   avg = hsi_conditional_average, inputvector[, limitfactor][, OK_FLAG=ok_flag]
;
; CALLS:
;   none
;
; INPUTS:
;   inputvector = nv element vector
;
; OPTIONAL INPUTS:
;   limitfactor = maximum allowable factor by which an input element can differ from the final average for it to contribute to that average.
;                   Default = 0.25
;
; OUTPUTS:
;   average value of input vector elements, excluding outliers
;
; OPTIONAL OUTPUTS:
;   none
;
; KEYWORDS:
;    ok_flag = nv-element vector with elements = 0 if corresponding element of the input vector element was rejected; element = 1 if ok.
;
; COMMON BLOCKS:
;   none
; SIDE EFFECTS:
;   none
;
; RESTRICTIONS:
;   For large quasi-gaussian distribution, outlier rejection could be obtained with hsi_outlier or resistant_mean.
;   This algorithm is better suited to smaller, non-gaussian distributions.
;
; MODIFICATION HISTORY:
;   8-Jan-2009  Initial version (ghurford@ssl.berkeley.edu)
;
;-
FUNCTION hsi_conditional_average, inputvector, limitfactor, OK_FLAG=ok_flag
;
DEFAULT, limitfactor, 0.25
nv          = N_ELEMENTS(inputvector)
ok_flag     = FLTARR(nv) + 1
indices     = FINDGEN(nv)
;
; To ensure algorithm stability, eliminate values one at a time.
FOR n=0, nv-1 DO BEGIN
    ok = indices[WHERE(ok_flag EQ 1, nok)]
    av = AVERAGE(inputvector[ok])
    IF nok LT 3 THEN BREAK                                                              ; exit loop if only 2 values left
    absdiff = ABS(inputvector[ok] - av)
    maxabsdiff = MAX(absdiff, imd)
    IF (maxabsdiff / ABS(av) GT limitfactor) THEN ok_flag[ok[imd]] = 0 ELSE BREAK       ; exit loop when no more outliers
ENDFOR
RETURN, av
END