##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
##MedFilter.py
##
##PURPOSE: 
##  Take a 1-d signal and return the median filtered version 
##
##INPUTS: 
##  Signal: float arr
##     The 1-d array of signal 
##  FilterSize: int
##     The full path to the slopefitted file corresponding to 
##  SigRej: float
##     The rejection value in terms of the standard deviation.  
##     Points SigRej*StdDev will not be included in the median
##  UseDJS: int
##     0 - Don't do sigma rejection
##     1 - Do sigma rejection  (NOTE: THIS IS VERY SLOW)
##  RejVal: float
##     A value that should not be included in the median calculation
##OUTPUTS:
##  MedSig: float arr
##    The median filtered signal
##  MeanSig: float arr
##    The mean filtered signal
##EXAMPLES:
##  from MedFilter import *
##  MedSig = MedFilter(Signal)
##
from Djs_Iterstat import *
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')

def MedFilter(Signal, FilterSize=8, SigRej=3.0, UseDJS=0, RejVal='null'):

  #Go through one example ramp and get the indices for the finite difference
  ZSize    = size(Signal)
  MedFilt  = zeros(ZSize, dtype=double)
  MeanFilt = zeros(ZSize, dtype=double)
  LoInds   = zeros(ZSize)
  HiInds   = zeros(ZSize)

  for ValInd in arange(ZSize):
    x0 = max((ValInd - FilterSize/2+1), 0)
    x1 = min((ValInd + FilterSize/2-1), ZSize-1)
    if x0 == 0       : x1 = ValInd+FilterSize/2
    if x1 > ZSize-1  : x0 = ValInd-FilterSize/2
    LoInds[ValInd]  = x0
    HiInds[ValInd]  = x1

  #Now go through the cube and make the estimates
  for k in arange(ZSize):
    x1     = HiInds[k]
    x0     = LoInds[k]
    SubSig = Signal[x0:x1] 

    #Determine whether or not to do DJS rejection
    if UseDJS == 1:
      Mean, Sig, Med, Mode, Mask =\
      Djs_Iterstat(SubSig , RejVal=RejVal, SigRej=SigRej, BinData=0)
    else:
      if RejVal != 'null' :
        GoodVals = (where(SubSig != RejVal))[0]
        if size(GoodVals) > 1 :
          Med  = median(SubSig[GoodVals])
          Mean = mean(SubSig[GoodVals])
        else:
          Med  = median(SubSig)
          Mean = Mean(SubSig)

    #Fill the array
    MedFilt[k]  = Med
    MeanFilt[k] = Mean 
 
  return MedFilt, MeanFilt