#===========================================================================
    # ;+
    # ; LOESS : Filter a signal H along X with a cutoff frequency CUT<br /><br />
    # ; 
    # ; Reference : Schlax, M. G., et D. B. Chelton (1992), Frequency Domain Diagnostics<br />
    # ;   for Linear Smoothers, Journal of the American Statistical Association, <br />
    # ;   87(420), 1070-1081. 
    # ; 
    # ; @param h {in}{required}{type:NUMERIC} Array of height anomalies
    # ; @param x {in}{required}{type:NUMERIC} Array of corresponding distances
    # ; @param cut {in}{required}{type:NUMERIC} Cutoff distance in units of X.
    # ; @param nan {in}{optional}{type:LOGICAL}{default:FALSE} Use NaNs<br />
    # ;          
    # ; @returns Hfilt (H filtered at cutoff CUT)
    # ; 
    # ; 
    # ; @uses exist
    # ; 
    # ; @author Renaud DUSSURGET, LEGOS/CTOH
    # ; @history Created Jan. 2010 by R.DUSSURGET from loess.c develloped by Mog2d team<br />
    # ;          May 2010 : Take NaNs into account for computation<br />
    # ;                     Checked the filters formulation (Schlax and Chelton, 1992) <br /><br />
    # ;
    # ;  LEGACY :<br />
    # ;  void dloess1d_irregular(int n, double l_c, double *h, double *x, double mask,double *lf,int *rstatus)<br /><br /> 
    # ;  
    # ;  ./varReduction.cpp:          dloess1d_irregular(n,cut,h,t,(double) satData.mask,lf,&status);<br />
    # ;  ./functions.extraction.cpp:  dloess1d_irregular(n,cut,h2,x,mask,lf,&status);<br /><br />
    # ;  
    # ;  rq: n-> nb points dans s�rie<br />
    # ;  l_c -> cut<br />
    # ;  *h -> sla<br />
    # ;  *x -> dimensional array<br />
    # ;  mask -> valeur du masque<br />
    # ;  *lf -> partie basse freq<br />
    # ;  *status -> retour<br />
    # ;-
    #===========================================================================

def loess(h,x,cut):
    """
    #===========================================================================
    # ;+
    # ; LOESS : Filter a signal H along X with a cutoff frequency CUT<br /><br />
    # ; 
    # ; Reference : Schlax, M. G., et D. B. Chelton (1992), Frequency Domain Diagnostics<br />
    # ;   for Linear Smoothers, Journal of the American Statistical Association, <br />
    # ;   87(420), 1070-1081. 
    # ; 
    # ; @param h {in}{required}{type:NUMERIC} Array of height anomalies
    # ; @param x {in}{required}{type:NUMERIC} Array of corresponding distances
    # ; @param cut {in}{required}{type:NUMERIC} Cutoff distance in units of X.
    # ; @param nan {in}{optional}{type:LOGICAL}{default:FALSE} Use NaNs<br />
    # ;          
    # ; @returns Hfilt (H filtered at cutoff CUT)
    # ; 
    # ; 
    # ; @uses exist
    # ; 
    # ; @author Renaud DUSSURGET, LEGOS/CTOH
    # ; @history Created Jan. 2010 by R.DUSSURGET from loess.c develloped by Mog2d team<br />
    # ;          May 2010 : Take NaNs into account for computation<br />
    # ;                     Checked the filters formulation (Schlax and Chelton, 1992) <br /><br />
    # ;
    # ;  LEGACY :<br />
    # ;  void dloess1d_irregular(int n, double l_c, double *h, double *x, double mask,double *lf,int *rstatus)<br /><br /> 
    # ;  
    # ;  ./varReduction.cpp:          dloess1d_irregular(n,cut,h,t,(double) satData.mask,lf,&status);<br />
    # ;  ./functions.extraction.cpp:  dloess1d_irregular(n,cut,h2,x,mask,lf,&status);<br /><br />
    # ;  
    # ;  rq: n-> nb points dans s�rie<br />
    # ;  l_c -> cut<br />
    # ;  *h -> sla<br />
    # ;  *x -> dimensional array<br />
    # ;  mask -> valeur du masque<br />
    # ;  *lf -> partie basse freq<br />
    # ;  *status -> retour<br />
    # ;-
    #===========================================================================
    """
          
    n=np.size(h)

    #If masked_array --> Get data and replace masked values by NaNs
    if isinstance(h,np.ma.masked_array): h=mask2NaN(h)
    else : h=np.ma.masked_array(h,mask=np.zeros(n,dtype=bool))

    #Start filtering if there are more than 1 point
    if n == 1 : lf = h[0]
    else : pass
    
    #Get the half-width of filter window
    l_c=cut/2.0          
    
    #Initiate weights and outputs
    w=np.zeros(n)
    lf = np.repeat(np.NaN,n)
    
    #Get valid points in serie
    flag = ~np.isnan(h)
    fcnt = flag.sum()
    fnt = np.arange(n).compress(flag)
    
    x = np.ma.masked_array(x,mask=(np.isnan(h) + h.mask) )
#        WHERE(FINITE(h),fcnt)

    #Vectorized form
    q = np.reshape(np.repeat(x,n),(n,n))
    q = (np.abs(q - q.transpose())/l_c).reshape(n*n)
    
    s = (1.0 - q*q*q)
    
#    outOfFilter_flag = q > 1.0
    outOfFilter_flag = np.greater(q,1)
    outCnt = np.nansum(outOfFilter_flag)
    outOfFilter = np.arange(n*n).compress(outOfFilter_flag)
    if (outCnt > 0) :
        s[outOfFilter]=0.0
        s.mask[outOfFilter]=True
    
    w=(s*s*s).reshape((n,n))  
    sumvar=np.nansum(w,1)

    #Compute current height (heights times weigths, divided by sum of weights)
    try :
        lf=np.nansum(w*np.repeat(h,n).reshape((n,n)),0)/sumvar
    except (ValueError):
        lf=np.nansum(w*np.repeat(h,n).reshape((n,n)),0)/sumvar #Retry

    return lf