;----------------------------------------------------------------------------
; $Id: scale_lyapunov.pro,v 1.13 2002/07/17 20:59:07 johnny Exp $
;+
; NAME:
; SCALE_LYAPUNOV
;
; PURPOSE:
; Calculate the scale-dependent Lyapunov exponent as defined by Hwang
; et al. (2000) given values of the time-dependent exponent (lambda).
;
; CATEGORY:
; Statistics.
;
; CALLING SEQUENCE:
; Result = SCALE_LYAPUNOV(in_lambda, in_npts, in_k, in_shell_bounds)
;
; INPUTS:
; in_lambda: Lambda for all in_shell_bounds intervals, where
; each row element i is the lambda in the interval:
; ( in_shell_bounds[i+1], in_shell_bounds[i] ),
; The NSH-1 lambda_all value is for the interval:
; ( 0, in_shell_bounds[NSH-1] ).
; Floating array dimensioned (NK,NSH). Unchanged by
; function.
;
; in_npts: Number of points used in each (k,shell) calculation.
; Each row is the value of the number of points for
; that shell, for all evolution times k. Long integer
; array dimensioned (NK,NSH). Unchanged by function.
;
; in_k: Evolution times (k) at which to calculate time-de-
; pendent exponent curves. Integer vector of NK ele-
; ments of equally spaced values of k, starting with
; 1. Unchanged by routine.
;
; in_shell_bounds: The shell upper boundaries for each in_lambda.
; The bounds for each i element in those output
; variables is the interval:
; ( in_shell_bounds(i+1), in_shell_bounds(i) ).
; The final element (NSH-1) corresponds to the
; interval:
; ( 0, in_shell_bounds(NSH-1) )
; which is a "shell" that is actually a "ball."
; Floating array dimensioned NSH. Unchanged by
; routine.
;
; KEYWORD PARAMETERS:
; AUTO_LIN_LIM: If keyword set true, specifies to use the automatic
; method of calculating "linear part" limits, and
; returns an array of the limits and number of points
; in the "linear part" via the Linear_Limits keyword.
; Either (but not both) keywords /Auto_Lin_Lim or
; /Manual_Lin_Lim must be true. See note below re-
; garding robustness of automatic method (which is
; very poor).
;
; IGNORE_L_MULT: If keyword is defined and set to the value of the
; delay time used to produce the lambda values, the
; function ignores k values that are integer mul-
; tiples of the delay time in calculating the Lya-
; punov exponents. This is to prevent artifactual
; correlations (see Cellucci et al. 1997).
;
; LINEAR_LIMITS: If Auto_Lin_Lim is true, this is an optional key-
; word, and is an output variable that creates/
; overwrites, returning a (3,NSH) integer array where
; the 1st col. is the lower (inclusive) limit k
; value, the 2nd col. is the upper (inclusive) limit
; k value, and the 3rd col. is the number of k values
; in the "linear part," for each shell. If a "linear
; part" doesn't exist, the lowe and upper limits for
; that shell are set to -1.
;
; If Manual_Lin_Lim is true, this is a required
; input keyword set to a (2,NSH) integer array,
; where the 1st col. is the lower limit k value,
; the 2nd col. is the upper limit k value, for each
; shell. All points between (and including) the
; lower and upper limits are considered in the
; "linear part" for each shell. If a "linear part"
; doesn't exist, set Linear_Limits to negative values,
; which will produce an NaN Lyapunov exponent.
;
; MANUAL_LIN_LIM: If keyword set true, specifies to use the manual
; method of calculating "linear part" limits, and
; reads an array of the limits of points in the
; "linear part" via the Linear_Limits keyword.
; Either (but not both) keywords /Auto_Lin_Lim or
; /Manual_Lin_Lim must be true.
;
; MIN_NPTS_IN_SHELL: If this is set to a variable, that variable is the
; value used as a test as to whether the value of
; lambda is reliable: if the number of points used
; to calculate lambda is less than this keyword's
; value, then lambda is set to NaN. Default value
; is 300. Variable unchanged by routine.
;
; OUTPUTS:
; Result: Scale-dependent Lyapunov exponent. Each element
; is the value of the scale-dependent Lyapunov
; exponent for a shell interval (see Shell_Bounds).
; If the Lyapunov exponent isn't well defined at
; a given scale, it is set to NaN. Floating array
; dimensioned NSH. Created by routine.
;
; FILE DEVICE I/O:
; None.
;
; COMMON BLOCKS:
; None.
;
; EXAMPLE:
; None.
;
; MODIFICATION HISTORY:
; - 31 May 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado.
; Email: air_jlin@yahoo.com. Passed minimally passably adequate tests.
; - 20 Jun 2002: Some small bug fixes for cases when there's no linear
; part defined.
; - 17 Jul 2002: Add ability to ignore multiples of delay time. Passed
; minimally passably adequate tests.
;
; NOTES:
; - Written for IDL 5.5.
; - References:
; * Cellucci, C. J., A. M. Albano, et al. (1997), "Detecting noise in a
; time series," Chaos, Vol. 7, No.3, pp. 414-422.
; * Gao, J. B. (1997), "Recognizing randomness in a time series," Physica
; D, Vol. 106, pp. 49-56
; * Hwang, S. K., J. B. Gao, and J. M. Liu (2000), "Noise-induced chaos
; in an optically injected semiconductor laser model," Phys. Rev. E,
; Vol. 61, No. 5, pp. 5162-5170
; - All keyword parameters are optional unless otherwise stated.
; - In the automated method of determining what part of the lambda curve
; is "linear," the key local variables are zerotest, nonzerotest,
; min_lin_k, and max_lin_k. I've tuned these values to work moderately
; well with one particular dataset; you may need to change them to fit
; your dataset.
; - NB: In comparing Auto_Lin_Lim with Manual_Lin_Lim, the automatic
; method is much poorer. I would recommend it never be used. I keep
; the keyword and code in the procedure, however, in case others wish
; to try to improve it in the future.
; - Negative Lyapunov exponents are not filtered out by the routine, but
; are returned as calculated.
; - No procedures called with _Extra keyword invoked.
; - No user-written procedures called.
;-
; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions,
; and contact information, see http://www.johnny-lin.com/lib.html.
;----------------------------------------------------------------------------
FUNCTION SCALE_LYAPUNOV, in_lambda, in_npts, in_k, in_shell_bounds $
, AUTO_LIN_LIM = auto_lin_lim $
, IGNORE_L_MULT = in_ignore_l_mult $
, LINEAR_LIMITS = io_linear_limits $
, MANUAL_LIN_LIM = manual_lin_lim $
, MIN_NPTS_IN_SHELL = in_min_npts_in_shell $
, _EXTRA = extra
; -------------------- Error Check and Parameter Setting --------------------
COMPILE_OPT IDL2
ON_ERROR, 0
lambda_all = in_lambda ;- protect input variables
npts_all = in_npts
k = in_k
shells_all = in_shell_bounds
NK = N_ELEMENTS(k) ;- set number of evolution times
NSH = N_ELEMENTS(shells_all) ;- set total number of shells
tmp = WHERE(SIZE(lambda_all, /Dimensions) ne [NK,NSH], count) ;- error chk.
if (count ne 0) then MESSAGE, 'error--bad in_lambda' ; key input
tmp = WHERE(SIZE(npts_all , /Dimensions) ne [NK,NSH], count) ; dims.
if (count ne 0) then MESSAGE, 'error--bad in_lambda'
if ((SIZE(k, /Type) ne 2) and $ ;- error check
(SIZE(k, /Type) ne 3)) then MESSAGE, 'error--bad k' ; k type
sum_check = KEYWORD_SET(auto_lin_lim) $ ;- check auto and manual
+ KEYWORD_SET(manual_lin_lim) ; keyword settings
if (sum_check ne 1) then $
MESSAGE, 'error--bad auto/man setting'
if (KEYWORD_SET(manual_lin_lim) eq 1) then begin ;- some checks and
if (N_ELEMENTS(io_linear_limits) eq 0) then $ ; settings for
MESSAGE, 'error--ll must be defined' ; /Manual_Lin_Lim case
linear_limits = io_linear_limits
sd_ll = SIZE(linear_limits)
tmp = WHERE(sd_ll[0:2] ne [2,2,NSH], count)
if (count ne 0) then MESSAGE, 'error--bad ll'
endif
if (N_ELEMENTS(in_min_npts_in_shell) gt 0) then begin ;- set minimum num.
min_npts_in_shell = in_min_npts_in_shell ; pts. a shell must
endif else begin ; have to have
min_npts_in_shell = 300 ; non-NaN lambda
endelse
if (N_ELEMENTS(in_ignore_l_mult) ne 0) then begin ;- set delay time that was
delay_time = in_ignore_l_mult ; used to generate the
endif ; lambda values
; ----------------- Find Largest Positive Lyapunov Exponent -----------------
;
; The slope of lambda (with evolution time as the abscissa) for the linear
; part ("small" evolution time) is the largest positive Lyapunov exponent.
; To calculate it we (1) determine where the "linear part" is, then (2)
; calculate the slope in that "linear part" by a "robust" least absolute
; deviation method (note that lambda must be sorted in x-ascending order
; before calculating the slope when using this method).
;
; There are two ways provided to do (1): (a) a manual method where the
; k values of the limits of the "linear part" are provided, or (b) an
; automatic method based on slopes and other criteria. Obviously, for
; (a) to work you have to visually inspect the lambda graphs before
; runtime. (a) will be more accurate than (b), however, because of the
; difficulty in creating an automatic method that filters out cases where
; there are spurious "linear parts" due to large oscillations in lambda.
;
; For (b), we determine the part that is linear by calculating the first
; and second derivatives at each point based on 3-point, Lagrangian
; interpolation. The linear part is that region where the second
; derivative is about zero and the first derivative is non-zero. To
; test these derivatives, we take the absolute value of the derivatives
; and normalize them by their maximum values and then test against
; zerotest and nonzerotest. We also restrict the linear part to values
; of k in the interval [min_lin_k, max_lin_k].
;
; The case statement in this section creates linear_pts and count_lp for
; the automatic and manual cases.
;
; If there are at least minfit number of points in the linear region, the
; Lyapunov exponent will be calculated. Key outputs of this section:
; lyapunov, linear_limits (in the /Auto_Lin_Lim case).
;
; Selected key variables:
; lyapunov Lyapunov exponent for each shell. Default is NaN,
; which is the value this is set to if a Lyapunov expo-
; nent isn't well-defined.
; count_lp Number of points in the "linear part." This value must
; be greater than minfit in order for the Lyapunov expo-
; nent to be calculated.
; max_lin_k Working variable for maximum and minimum bounds for
; min_lin_k "linear part." In the automatic case, these are just
; the outer brackets; most of the time the actual bounds
; defining the "linear part" are within these max./min.
; For the manual case, these are set to the bounds that
; were input into this routine.
; minfit Min. num. pts. needed in "linear range" to calculate
; slope for Lyapunov exponent.
; nonzerotest The test value for whether a normalized absolute value
; is non-zero; values more than nonzerotest are considered
; non-zero.
; zerotest The test value for whether a normalized absolute value
; is zero; values less than zerotest are considered zero.
if (KEYWORD_SET(auto_lin_lim) eq 1) then begin ;- create linear_limits
linear_limits = REPLICATE(-1,3,NSH) ; for /Auto_Lin_Lim case
MESSAGE, 'warning--are you sure you want ' $
+ 'to use the auto method? comment ' $
+ 'this line out if so.'
endif
lyapunov = REPLICATE(!VALUES.F_NAN, NSH) ;- est. Lyap. expon. array
minfit = 3 ;- min. pts. needed in "lin-
; ear range" for calcs.
for ishell=0,NSH-1 do begin ;-> cycle through shells (begin)
lambda = REFORM(lambda_all[*,ishell]) ;- lambda for this ishell
case 1 of ;-> choose between auto and manual settings (begin)
KEYWORD_SET(auto_lin_lim): begin
zerotest = 0.10 ;- pts. less than this value are "0"
nonzerotest = 0.20 ;- pts. more than this value are "non-0"
min_lin_k = 1 ;- min. k value for the linear part
max_lin_k = 70 ;- max. k value for the linear part
slope = DERIV(k,lambda) ;- calc. abs. value
absslope = ABS( slope ) ; of 1st and 2nd
absderiv2 = ABS( DERIV(k,slope) ) ; deriv. and the
absslope_norm = absslope/MAX(absslope) ; linear part pts.
absderiv2_norm = absderiv2/MAX(absderiv2)
linear_pts = WHERE( (absslope_norm gt nonzerotest) and $
(absderiv2_norm lt zerotest) and $
(k ge min_lin_k) and $
(k le max_lin_k), count_lp )
if (count_lp gt 0) then begin ;- fill array
linear_limits[0,ishell] = k[MIN(linear_pts)] ; linear_limits
linear_limits[1,ishell] = k[MAX(linear_pts)]
endif
linear_limits[2,ishell] = count_lp
end
KEYWORD_SET(manual_lin_lim): begin
min_lin_k = linear_limits[0,ishell] ;- min. k, linear part
max_lin_k = linear_limits[1,ishell] ;- max. k, linear part
linear_pts = WHERE( (k ge min_lin_k) and $
(k le max_lin_k), count_lp )
end
else: MESSAGE, 'error--bad auto/man settings'
endcase ;-> choose between auto and manual settings (end)
if ( (N_ELEMENTS(in_ignore_l_mult) ne 0) and $ ;- ignore k's that
(count_lp ge minfit) ) then begin ; are multiples
not_L_mult_pts = $ ; of delay_time
WHERE( (k[linear_pts] mod delay_time) ne 0 $
, count )
if (count gt 0) then linear_pts=linear_pts[not_L_mult_pts]
count_lp = N_ELEMENTS(linear_pts)
endif
if (count_lp ge minfit) then begin ;- if there's enough pts.,
kf = k[linear_pts] ; calc. slope of lambda
lf = lambda[linear_pts] ; in the "linear part",
kf = kf[SORT(kf)] ; for pts. sorted in x-
lf = lf[SORT(kf)] ; ascending order
lyapunov[ishell] = (LADFIT(kf, lf))[1]
endif
endfor ;-> cycle through shells (end)
; ------------------------- Check for Npts Too Small ------------------------
;
; The number of points in each shell used to calculate lambda is tested
; against Min_Npts_In_Shell to see whether one can have confidence in the
; Lyapunov value derived from lambda (if the number of points is less than
; min_npts_in_shell, lyapunov is set to NaN).
;
; If npts_all is not constant for all k's for a given shell, a warning
; message is printed, and the value of npts_all used in the test against
; min_npts_in_shell is the minimum of npts_all in that shell. Those
; minimum values are stored in min_npts_all.
min_npts_all = LONARR(NSH) ;- minimum of npts_all in each shell
for ishell=0,NSH-1 do begin ;- check npts is identical
tmpnpts = npts_all[*,ishell] ; for all k's for a given
if (MAX(tmpnpts) ne MIN(tmpnpts)) then $ ; shell
PRINT, 'warning--npts not same'
min_npts_all[ishell] = MIN(tmpnpts)
endfor
bad_pts = WHERE( min_npts_all lt $ ;- find unreliable
min_npts_in_shell, count ) ; points and set
if (count gt 0) then lyapunov[bad_pts] = !VALUES.F_NAN ; to NaN
; ------------------------------ Prepare Output -----------------------------
;
; Return Lyapunov exponent, shell bounds, and limits of the linear part.
if (KEYWORD_SET(auto_lin_lim) eq 1) then begin ;- set limits of linear
io_linear_limits = linear_limits ; part for /Auto_Lin_Lim
endif ; (optional)
Result = TEMPORARY(lyapunov) ;- return scale-dependent
; Lyapunov exponent
RETURN, Result
END ;=== end of function ===
; ========== end file ==========