Index: trunk/SRC/ToBeReviewed/STATISTICS/a_timecorrelate.pro
===================================================================
--- trunk/SRC/ToBeReviewed/STATISTICS/a_timecorrelate.pro	(revision 334)
+++ trunk/SRC/ToBeReviewed/STATISTICS/a_timecorrelate.pro	(revision 335)
@@ -16,11 +16,11 @@
 ; @keyword ZERO2NAN
 ;
+; @keyword NAN
+;
 ; @keyword DOUBLE
 ; If set to a non-zero value, computations are done in
 ; double precision arithmetic.
 ;
-; @examples
-;
-; @history
+; @hidden
 ;
 ; @version
@@ -28,33 +28,37 @@
 ;
 ;-
-FUNCTION timeauto_cov, x, m, nt, DOUBLE=double, ZERO2NAN=zero2nan
+FUNCTION timeauto_cov, x, m, nt, DOUBLE=double, ZERO2NAN=zero2nan, NAN =  nan
+;Sample autocovariance function
 ;
   compile_opt idl2, strictarrsubs
 ;
-;Sample autocovariance function
-   TimeDim = size(X, /n_dimensions)
-   Xmean = TOTAL(X, TimeDim, Double = Double) / nT
-   if double then one = 1.0d ELSE one = 1.0
-   Xmean = Xmean[*]#replicate(one, nT - M)
-;
-;
-   case TimeDim of
-      1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (X[M:nT - 1L] - Xmean), $
-                    TimeDim, Double = Double)
-      2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * $
-                    (X[*, M:nT - 1L] - Xmean) $
-                    , TimeDim, Double = Double)
-      3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * $
-                    (X[*, *, M:nT - 1L] - Xmean) $
-                    , TimeDim, Double = Double)
-      4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * $
-                    (X[*, *, *, M:nT - 1L] - Xmean) $
-                    , TimeDim, Double = Double)
-   ENDCASE
-   if keyword_set(zero2nan) then begin
-      zero = where(res EQ 0)
-      if zero[0] NE -1 then res[zero] = !values.f_nan
-   endif
-   RETURN, res
+  IF NAN AND M GE 1 THEN $
+   STOP, 'Warning : lagged autocorrelation is not possible at the moment for time-series with NaN !!!'
+
+  TimeDim = size(X, /n_dimensions)
+  Xmean = NAN ? TOTAL(X, TimeDim, Double = Double,  NAN =  nan) / TOTAL(FINITE(X), TimeDim) : $
+   TOTAL(X, TimeDim, Double = Double) / nT                                                                                
+  one = double ? 1.0d : 1.0
+  Xmean = Xmean[*]#replicate(one, nT - M)
+
+; Time-series with NaN : only for Lag = 0
+  case TimeDim of
+     1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (X[M:nT - 1L] - Xmean), $
+                   TimeDim, Double = Double, NAN =  nan)
+     2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * $
+                   (X[*, M:nT - 1L] - Xmean[*]) $
+                   , TimeDim, Double = Double, NAN =  nan)
+     3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean[*, *]) * $
+                   (X[*, *, M:nT - 1L] - Xmean) $
+                   , TimeDim, Double = Double, NAN =  nan)
+     4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * $
+                   (X[*, *, *, M:nT - 1L] - Xmean) $
+                   , TimeDim, Double = Double,  NAN =  nan)
+  ENDCASE
+  if keyword_set(zero2nan) then begin
+     zero = where(res EQ 0)
+     if zero[0] NE -1 then res[zero] = !values.f_nan
+  endif
+  RETURN, res
 
 END
@@ -78,6 +82,6 @@
 ;
 ; @param LAG {in}{required}{type=scalar or vector}
-; A scalar or n-elements vector, in the interval [-(n-2),(n-2)],
-; of type integer that specifies the absolute distance(s) between
+; A scalar or n-element vector, in the interval [-(n-2), (n-2)],
+; of type integer that specifies the absolute distance(s) between 
 ; indexed elements of X.
 ;
@@ -86,4 +90,8 @@
 ; is computed.
 ;
+; @keyword NVAL
+; A named variable that, on exit, contains the number of valid
+; observations (not NAN)
+;
 ; @keyword DOUBLE
 ; If set to a non-zero value, computations are done in
@@ -91,5 +99,5 @@
 ;
 ; @examples
-; Define an n-elements sample population.
+; Define an n-element sample population.
 ; IDL> x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
 ;
@@ -113,5 +121,5 @@
 ;
 ;-
-FUNCTION a_timecorrelate, x, lag, COVARIANCE=covariance, DOUBLE=double
+FUNCTION a_timecorrelate, x, lag, COVARIANCE=covariance, DOUBLE=double, NVAL =  nval
 ;
   compile_opt idl2, strictarrsubs
@@ -126,13 +134,22 @@
    XNDim = SIZE(X, /n_dimensions)
    nT = XDim[XNDim-1]
-                                ;Check length.
-   if nT lt 2 then $
-    ras= report("Time axis of X array must contain 2 or more elements.")
+
+; Keyword NAN activated if needed for TimeAuto_Cov function
+; Keyword NVAL not compulsory.
+
+   NAN = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ? 1 : 0
+;We can retrieve the matrix of real lenghts of time-series   
+   nTreal = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ?  TOTAL(FINITE(X),  XNDim) : nT 
+
+   IF ARG_PRESENT(NVAL) THEN nval = nTreal
+
+;Check length.
+   IF (WHERE(nTreal LE 1))[0] NE -1 THEN $
+    MESSAGE,  "Matrix of length of time-series must contain 2 or more elements"   
 
 ;If the DOUBLE keyword is not set then the internal precision and
 ;result are identical to the type of input.
-   if N_ELEMENTS(Double) eq 0 then $
-    Double = (SIZE(X, /type) eq 5)
-
+   type = SIZE(X, /TYPE)
+   useDouble = (N_ELEMENTS(Double) eq 1) ? KEYWORD_SET(Double) : (type eq 5)
 
    if n_elements(lag) EQ 0 then lag = 0
@@ -141,37 +158,36 @@
    if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
 
+; Type of outputs according to the type of data input
+
    case XNDim of
-      1:if Double eq 0 then Auto = FLTARR(nLag) else Auto = DBLARR(nLag)
-      2:if Double eq 0 then Auto = FLTARR(XDim[0], nLag) else Auto = DBLARR(XDim[0], nLag)
-      3:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], nLag) $
-      else Auto = DBLARR(XDim[0], XDim[1], nLag)
-      4:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], XDim[2], nLag) $
-      else Auto = DBLARR(XDim[0], XDim[1], XDim[2], nLag)
+      1: Auto = useDouble ? DBLARR(nLag) : FLTARR(nLag)
+      2: Auto = useDouble ? DBLARR(XDim[0], nLag) : FLTARR(XDim[0], nLag)
+      3: Auto = useDouble ? DBLARR(XDim[0], XDim[1], nLag) : FLTARR(XDim[0], XDim[1], nLag)
+      4: Auto = useDouble ? DBLARR(XDim[0], XDim[1], XDim[2], nLag) : FLTARR(XDim[0], XDim[1], XDim[2], nLag)
    endcase
 
-   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Autocorrelation.
-      for k = 0, nLag-1 do $
-       case XNDim of
-         1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
-          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
-         2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
-          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
-         3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
-          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
-         4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
-          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
-      endcase
-   endif else begin             ;Compute Autocovariance.
-      for k = 0, nLag-1 do $
-       case XNDim of
-         1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
-         2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
-         3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
-         4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
-      endcase
-   endelse
+; Compute lagged autocorrelation or autocovariance (no NaN)
+   FOR k = 0, nLag-1 DO BEGIN
+      case XNDim of
+         1: BEGIN
+            Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $
+             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble,  /zero2nan, NAN = nan) )
+         END
+         2: BEGIN
+            Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $
+             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) )   
+         END
+         3: BEGIN
+            Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $
+             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) )
+         END
+         4: BEGIN
+            Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $
+             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) )
+         END
+      ENDCASE
+   ENDFOR
 
-   if Double eq 0 then RETURN, FLOAT(Auto) else $
-    RETURN, Auto
+   return, useDouble ? Auto : FLOAT(Auto)
 
 END
Index: trunk/SRC/ToBeReviewed/STATISTICS/skewness_4d.pro
===================================================================
--- trunk/SRC/ToBeReviewed/STATISTICS/skewness_4d.pro	(revision 335)
+++ trunk/SRC/ToBeReviewed/STATISTICS/skewness_4d.pro	(revision 335)
@@ -0,0 +1,137 @@
+;+
+;
+; @file_comments
+; 
+; @categories
+; Statistics
+;
+; @param X {in}{required}{type=array}
+; An Array which last dimension is the time dimension so
+; size n.
+;
+; @param NT
+;
+; @keyword DOUBLE
+; If set to a non-zero value, computations are done in
+; double precision arithmetic.
+;
+; @hidden
+;
+; @version
+; $Id$
+;
+;-
+FUNCTION Skewness_Num, X, nT, Double = Double,  NAN =  nan
+; Compute the numerator of the skewness expression
+;
+
+  compile_opt idl2, strictarrsubs
+;  
+  TimeDim = size(X, /n_dimensions)
+  Xmean = NAN ? TOTAL(X, TimeDim, Double = Double,  NAN =  nan) / TOTAL(FINITE(X), TimeDim) : $
+   TOTAL(X, TimeDim, Double = Double) / nT    
+  one = double ? 1.0d : 1.0
+  Xmean = Xmean[*]#replicate(one, nT)
+  res = TOTAL( (X-Xmean)^3, TimeDim, Double = Double, NAN =  nan)
+
+  RETURN, res
+
+END
+;+
+; @file_comments
+; Same function as SKEWNESS but accept array (until 4
+; dimension) for input and perform  the skewness 
+; along the time dimension which must be the last
+; one of the input array.
+;
+; @categories
+; Statistics
+;
+; @param X {in}{required}{type=array}
+; An Array which last dimension is the time dimension so
+; size n.
+;
+; @keyword DOUBLE
+; If set to a non-zero value, computations are done in
+; double precision arithmetic.
+;
+; @keyword NVAL
+; A named variable that, on exit, contains the number of valid
+; observations (not NAN)
+;
+; @examples
+;
+; @history
+; May 2007 Michel Kolasinski (michel.kolasinski@locean-ipsl.upmc.fr)
+;
+; Based on the  a_timecorrelate procedure of IDL
+; INTRODUCTION TO STATISTICAL TIME SERIES
+; Wayne A. Fuller
+; ISBN 0-471-28715-6
+;
+; @version
+; $Id$
+;
+;-
+
+FUNCTION skewness_4d, X, DOUBLE = Double,  NVAL =  nval
+;
+   compile_opt idl2, strictarrsubs
+;
+; Compute the skewness from 1d to 4d vectors
+
+   ON_ERROR, 2
+
+   XDim = SIZE(X, /dimensions)
+   XNDim = SIZE(X, /n_dimensions)
+   nT = XDim[XNDim-1]
+
+; Keyword NAN activated if needed
+; Keyword NVAL not compulsory.
+
+   NAN = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ? 1 : 0
+;We can retrieve the matrix of real lenghts of time-series   
+   nTreal = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ?  TOTAL(FINITE(X),  XNDim) : nT 
+
+   IF ARG_PRESENT(NVAL) THEN nval = nTreal
+
+; Check length.
+   IF (WHERE(nTreal LE 1))[0] NE -1 THEN $
+    MESSAGE,  "Matrix of length of time-series must contain 2 or more elements"   
+
+; If the DOUBLE keyword is not set then the internal precision and
+; result are identical to the type of input.
+   type = SIZE(X, /TYPE)
+   useDouble = (N_ELEMENTS(Double) eq 1) ? KEYWORD_SET(Double) : $
+    (type eq 5)
+
+; Type of outputs according to the type of data input
+   case XNDim of
+      1: Skew = useDouble ? DBLARR(1) : FLTARR(1)
+      2: Skew = useDouble ? DBLARR(XDim[0]) : FLTARR(XDim[0])
+      3: Skew = useDouble ? DBLARR(XDim[0], XDim[1]) : FLTARR(XDim[0], XDim[1])
+      4: Skew = useDouble ? DBLARR(XDim[0], XDim[1], XDim[2]) : FLTARR(XDim[0], XDim[1], XDim[2])
+   endcase
+; Compute standard deviation
+; nTreal might be a matrix
+   std = a_timecorrelate(X, 0, /covariance)
+   std = sqrt(std)
+   zero = where(std EQ 0)
+
+   if zero[0] NE -1 then STOP,  $
+    'Cannot compute skewness since there are zeros in the matrix of standard deviations !'
+
+; Problem with high masked values (x^3 makes NaN when x is high)
+; Threshold put on the values of the tab
+   idx_std = WHERE (std GT 1.0e+10)
+   X = X < 1.0e+10
+   std = std < 1.0e+10
+
+; Compute skewness
+   Skew = Skewness_Num(X, nT, Double = useDouble, NAN = nan) / (nTreal*std^3)
+   IF idx_std[0] NE -1 THEN Skew[idx_std] = valmask
+
+   return, useDouble ? Skew : FLOAT(Skew)
+
+END
+
