;---------------------------------------------------------------------------- ; $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 ==========