Changeset 335 for trunk


Ignore:
Timestamp:
03/07/08 16:09:25 (16 years ago)
Author:
smasson
Message:

add Michel developments

Location:
trunk/SRC/ToBeReviewed/STATISTICS
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/ToBeReviewed/STATISTICS/a_timecorrelate.pro

    r327 r335  
    1616; @keyword ZERO2NAN 
    1717; 
     18; @keyword NAN 
     19; 
    1820; @keyword DOUBLE 
    1921; If set to a non-zero value, computations are done in 
    2022; double precision arithmetic. 
    2123; 
    22 ; @examples 
    23 ; 
    24 ; @history 
     24; @hidden 
    2525; 
    2626; @version 
    … …  
    2828; 
    2929;- 
    30 FUNCTION timeauto_cov, x, m, nt, DOUBLE=double, ZERO2NAN=zero2nan 
     30FUNCTION timeauto_cov, x, m, nt, DOUBLE=double, ZERO2NAN=zero2nan, NAN =  nan 
     31;Sample autocovariance function 
    3132; 
    3233  compile_opt idl2, strictarrsubs 
    3334; 
    34 ;Sample autocovariance function 
    35    TimeDim = size(X, /n_dimensions) 
    36    Xmean = TOTAL(X, TimeDim, Double = Double) / nT 
    37    if double then one = 1.0d ELSE one = 1.0 
    38    Xmean = Xmean[*]#replicate(one, nT - M) 
    39 ; 
    40 ; 
    41    case TimeDim of 
    42       1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (X[M:nT - 1L] - Xmean), $ 
    43                     TimeDim, Double = Double) 
    44       2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * $ 
    45                     (X[*, M:nT - 1L] - Xmean) $ 
    46                     , TimeDim, Double = Double) 
    47       3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * $ 
    48                     (X[*, *, M:nT - 1L] - Xmean) $ 
    49                     , TimeDim, Double = Double) 
    50       4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * $ 
    51                     (X[*, *, *, M:nT - 1L] - Xmean) $ 
    52                     , TimeDim, Double = Double) 
    53    ENDCASE 
    54    if keyword_set(zero2nan) then begin 
    55       zero = where(res EQ 0) 
    56       if zero[0] NE -1 then res[zero] = !values.f_nan 
    57    endif 
    58    RETURN, res 
     35  IF NAN AND M GE 1 THEN $ 
     36   STOP, 'Warning : lagged autocorrelation is not possible at the moment for time-series with NaN !!!' 
     37 
     38  TimeDim = size(X, /n_dimensions) 
     39  Xmean = NAN ? TOTAL(X, TimeDim, Double = Double,  NAN =  nan) / TOTAL(FINITE(X), TimeDim) : $ 
     40   TOTAL(X, TimeDim, Double = Double) / nT                                                                                 
     41  one = double ? 1.0d : 1.0 
     42  Xmean = Xmean[*]#replicate(one, nT - M) 
     43 
     44; Time-series with NaN : only for Lag = 0 
     45  case TimeDim of 
     46     1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (X[M:nT - 1L] - Xmean), $ 
     47                   TimeDim, Double = Double, NAN =  nan) 
     48     2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * $ 
     49                   (X[*, M:nT - 1L] - Xmean[*]) $ 
     50                   , TimeDim, Double = Double, NAN =  nan) 
     51     3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean[*, *]) * $ 
     52                   (X[*, *, M:nT - 1L] - Xmean) $ 
     53                   , TimeDim, Double = Double, NAN =  nan) 
     54     4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * $ 
     55                   (X[*, *, *, M:nT - 1L] - Xmean) $ 
     56                   , TimeDim, Double = Double,  NAN =  nan) 
     57  ENDCASE 
     58  if keyword_set(zero2nan) then begin 
     59     zero = where(res EQ 0) 
     60     if zero[0] NE -1 then res[zero] = !values.f_nan 
     61  endif 
     62  RETURN, res 
    5963 
    6064END 
    … …  
    7882; 
    7983; @param LAG {in}{required}{type=scalar or vector} 
    80 ; A scalar or n-elements vector, in the interval [-(n-2),(n-2)], 
    81 ; of type integer that specifies the absolute distance(s) between 
     84; A scalar or n-element vector, in the interval [-(n-2), (n-2)], 
     85; of type integer that specifies the absolute distance(s) between  
    8286; indexed elements of X. 
    8387; 
    … …  
    8690; is computed. 
    8791; 
     92; @keyword NVAL 
     93; A named variable that, on exit, contains the number of valid 
     94; observations (not NAN) 
     95; 
    8896; @keyword DOUBLE 
    8997; If set to a non-zero value, computations are done in 
    … …  
    9199; 
    92100; @examples 
    93 ; Define an n-elements sample population. 
     101; Define an n-element sample population. 
    94102; IDL> x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57] 
    95103; 
    … …  
    113121; 
    114122;- 
    115 FUNCTION a_timecorrelate, x, lag, COVARIANCE=covariance, DOUBLE=double 
     123FUNCTION a_timecorrelate, x, lag, COVARIANCE=covariance, DOUBLE=double, NVAL =  nval 
    116124; 
    117125  compile_opt idl2, strictarrsubs 
    … …  
    126134   XNDim = SIZE(X, /n_dimensions) 
    127135   nT = XDim[XNDim-1] 
    128                                 ;Check length. 
    129    if nT lt 2 then $ 
    130     ras= report("Time axis of X array must contain 2 or more elements.") 
     136 
     137; Keyword NAN activated if needed for TimeAuto_Cov function 
     138; Keyword NVAL not compulsory. 
     139 
     140   NAN = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ? 1 : 0 
     141;We can retrieve the matrix of real lenghts of time-series    
     142   nTreal = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ?  TOTAL(FINITE(X),  XNDim) : nT  
     143 
     144   IF ARG_PRESENT(NVAL) THEN nval = nTreal 
     145 
     146;Check length. 
     147   IF (WHERE(nTreal LE 1))[0] NE -1 THEN $ 
     148    MESSAGE,  "Matrix of length of time-series must contain 2 or more elements"    
    131149 
    132150;If the DOUBLE keyword is not set then the internal precision and 
    133151;result are identical to the type of input. 
    134    if N_ELEMENTS(Double) eq 0 then $ 
    135     Double = (SIZE(X, /type) eq 5) 
    136  
     152   type = SIZE(X, /TYPE) 
     153   useDouble = (N_ELEMENTS(Double) eq 1) ? KEYWORD_SET(Double) : (type eq 5) 
    137154 
    138155   if n_elements(lag) EQ 0 then lag = 0 
    … …  
    141158   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector. 
    142159 
     160; Type of outputs according to the type of data input 
     161 
    143162   case XNDim of 
    144       1:if Double eq 0 then Auto = FLTARR(nLag) else Auto = DBLARR(nLag) 
    145       2:if Double eq 0 then Auto = FLTARR(XDim[0], nLag) else Auto = DBLARR(XDim[0], nLag) 
    146       3:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], nLag) $ 
    147       else Auto = DBLARR(XDim[0], XDim[1], nLag) 
    148       4:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], XDim[2], nLag) $ 
    149       else Auto = DBLARR(XDim[0], XDim[1], XDim[2], nLag) 
     163      1: Auto = useDouble ? DBLARR(nLag) : FLTARR(nLag) 
     164      2: Auto = useDouble ? DBLARR(XDim[0], nLag) : FLTARR(XDim[0], nLag) 
     165      3: Auto = useDouble ? DBLARR(XDim[0], XDim[1], nLag) : FLTARR(XDim[0], XDim[1], nLag) 
     166      4: Auto = useDouble ? DBLARR(XDim[0], XDim[1], XDim[2], nLag) : FLTARR(XDim[0], XDim[1], XDim[2], nLag) 
    150167   endcase 
    151168 
    152    if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Autocorrelation. 
    153       for k = 0, nLag-1 do $ 
    154        case XNDim of 
    155          1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $ 
    156           TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan) 
    157          2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $ 
    158           TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan) 
    159          3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $ 
    160           TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan) 
    161          4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $ 
    162           TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan) 
    163       endcase 
    164    endif else begin             ;Compute Autocovariance. 
    165       for k = 0, nLag-1 do $ 
    166        case XNDim of 
    167          1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT 
    168          2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT 
    169          3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT 
    170          4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT 
    171       endcase 
    172    endelse 
     169; Compute lagged autocorrelation or autocovariance (no NaN) 
     170   FOR k = 0, nLag-1 DO BEGIN 
     171      case XNDim of 
     172         1: BEGIN 
     173            Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $ 
     174             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble,  /zero2nan, NAN = nan) ) 
     175         END 
     176         2: BEGIN 
     177            Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $ 
     178             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) )    
     179         END 
     180         3: BEGIN 
     181            Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $ 
     182             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) ) 
     183         END 
     184         4: BEGIN 
     185            Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = useDouble, NAN = nan) / $ 
     186             ( KEYWORD_SET(Covariance) ?  nTreal : TimeAuto_Cov(X, 0L, nT, Double = useDouble, /zero2nan, NAN = nan) ) 
     187         END 
     188      ENDCASE 
     189   ENDFOR 
    173190 
    174    if Double eq 0 then RETURN, FLOAT(Auto) else $ 
    175     RETURN, Auto 
     191   return, useDouble ? Auto : FLOAT(Auto) 
    176192 
    177193END 
Note: See TracChangeset for help on using the changeset viewer.