!============================================================================ ! General Documentation Section ! Name: ! LAMBDA_LOG_DISP ! ! Purpose: ! Calculate the time-dependent exponent (lambda) and logarithmic dis- ! placement curves, as defined by Gao (1997). ! ! Calling Sequence: ! call LAMBDA_LOG_DISP( x0, k, m, L, ie_value & ! , tan_motion_lim & ! , normalize_input & ! , out_lambda, out_log_disp & ! , out_npts, out_shells ) ! ! This subroutine must have an explicit interface where called. ! ! Example: ! None. ! ! File I/O: ! None. ! ! Entry and Exit States: ! N/A. ! ! RCS Information: ! $Id: lambda_log_disp.f90,v 1.1 2003/04/04 04:39:27 jlin Exp jlin $ ! ! Revision History: ! - 1993: Original version by Jianbo Gao, UCLA Dept. of Electrical Engin- ! eering. ! - 22 May 2002: Fortran 90 commenting and formatting by Johnny Lin, ! CIRES/University of Colorado, Boulder. Email: air_jlin@yahoo.com. ! Passed a few minimal tests. ! ! Notes: ! - References: ! * Gao, J. B. (1997), "Recognizing randomness in a time series," ! Physica D, Vol. 106, pp. 49-56. ! * Gao J. B., and Zheng Z. M. (1993), "Local exponential divergence ! plot and optimal embedding of a chaotic time-series," Phys. Lett. A, ! Vol. 181, No. 2, pp. 153-158. ! * Gao, J. and Z. Zheng (1994), "Direct dynamical test for deterministic ! chaos and optimal embedding of a chaotic time series," Phys. Rev. E, ! Vol. 49, No. 5, pp. 3807-3814. ! * Theiler, J. (1990), "Estimating fractal dimension," J. Opt. Soc. Am. ! A, Vol. 7, No. 6, pp. 1055-1073. ! - See http://www.johnny-lin.com/fort_refs/lambda_log_disp.notes.html for ! important notes regarding parameter choice. ! - In computations, normalized time series (i.e. in interval [0,1]) is ! assumed. If the input is not normalized, set parameter normalize_input ! to 1 to make timeseries normalized before beginning computations. ! - Prescribed rulers (squared): Given an array rl(9), of base 2 and ! exponent (-i) {i=ie,ie+1,...}, to define shells. Note that to properly ! define the ensemble averages at least several hundreds pairs of points ! are needed. See out_npts for the numbers used in averaging. ! - The number of points actually used in calculation (local variable nn) ! is not the same as the number of points in the input timeseries, ! because one does not necessarily need all those points to adequately ! span the phase space. Thus, "nn" is set to a value based on input ! timeseries length, embedding dimension, etc. or 8000, whichever is ! smaller. If this is unacceptable, goto the "Main Parameters Setting" ! sub-section and change "nn" accordingly. ! - Note: This routine does not use "implicit none," because when it ! was first written, there were a large number of variables that were ! implicitly declared. ! - No subprograms required by this routine besides those included in ! the Fortran 90 standard. ! ! Copyright (c) 1993-2002 by Jianbo Gao and Johnny Lin. For licensing ! and contact information see http://www.johnny-lin.com/flib.html. !---------------------------------------------------------------------------- ! Variable Declaration and Description Section ! Overall Routine Declaration: subroutine lambda_log_disp( x0, k, m, L, ie_value & , tan_motion_lim & , normalize_input & , out_lambda, out_log_disp & , out_npts, out_shells ) ! Include Files: ! None. ! Modules Used and Module Variables Changed by Subroutine: ! None. ! Input Arguments: real, intent(in) & !- Input scalar timeseries. Real vector. :: x0(:) integer, intent(in) & !- Evolution times at which to calculate :: k(:) ! time-dependent exponent. Integer ! vector of equally spaced values of k, ! starting with 1. integer, intent(in) & !- Embedding dimension. Scalar integer. :: m integer, intent(in) & !- Delay time. Scalar integer. :: L integer, intent(in) & !- Ruler initialization (to give rl(1)). :: ie_value ! This variable defines shells used in ! calculations and output in out_shells; ! each shell is defined as the interval: ! ( 2**[-(i+1)/2], 2**[-i/2] ) ! where i is an integer from ie_value to ! ie_value+8. integer, intent(in) & !- Value of the criteria to remove tangen- :: tan_motion_lim ! tial motion effects (i.e. ABS(i-j) must ! be greater than this value to be consid- ! ered for calculations). If parameter ! is set to < 0, the subroutine uses a ! value of 10 for the tangential effects ! removal test. integer, intent(in) & !- If value is 1, before using the input :: normalize_input ! time0 series x0, the subroutine normal- ! izes it by: (x0 - min)/(max - min), ! where min and max are the maximum and ! minimum of x0, respectively (i.e. normal- ! ize to unit interval [0,1]). ! Output Arguments: real, intent(out) & !- Each row is the time-dependent expon- :: out_lambda(SIZE(k),9) ! ent (lambda) at each evolution time k, ! for shells given by the shell interval ! (out_shells(i+1), out_shells(i)). real, intent(out) & !- Each row is the logarithmic displace- :: out_log_disp(SIZE(k),9) ! ment at each evolution time k, for ! shells given by the shell interval ! (out_shells(i+1), out_shells(i)). integer, intent(out) & !- Number of point pairs used for the en- :: out_npts(SIZE(k),9) ! semble averaging in each shell at each ! evolution time k, for shells given ! by the shell interval (out_shells(i+1), ! out_shells(i)). real, intent(out) & !- Shell upper boundaries for each row :: out_shells(9) ! in out_lambda and out_log_disp. The ! bounds for each i row in those output ! variables is the interval: ! ( out_shells(i+1), out_shells(i) ). ! The final row (i=9) corresponds to ! the interval: ! ( 0, out_shells(9) ) ! which is a "shell" that is actually a ! "ball." ! Local Variables: real, dimension(SIZE(k)+1,50) :: a real, dimension(SIZE(k)+1) :: d real, dimension(SIZE(x0)) :: xy real, dimension(SIZE(x0)) :: x !- normalized ver. of x0 real :: r(SIZE(k)+1) real :: ssz(SIZE(k),9) real :: rl(9) integer :: mapped(SIZE(k)) integer :: nsz(SIZE(k),9) integer :: nxp !- num pts. in x0 time series integer :: map !- num pts. in evol. time vector integer :: nxp_p, nxp1, map1 real :: cmax, cmin, ci integer :: tan_motion_lim0 integer :: nn !- num pts. actually used in calcs. integer :: ie, in, md, mt integer :: mstep, mmap integer :: i !============================================================================ ! ------------------------ Initial Parameter Setting ------------------------ if (tan_motion_lim .ge. 0) then !- set test to remove tangential motion tan_motion_lim0 = tan_motion_lim ! effects else tan_motion_lim0 = 10 endif nxp = SIZE(x0) !- set some array bounds etc. nxp_p = SIZE(x0) nxp1 = SIZE(x0) map = SIZE(k) map1 = map + 1 x = x0 !- set x to x0 (normalize below, if chosen) ! ----------------------- Normalize Input (If Chosen) ----------------------- if (normalize_input .eq. 1) then cmax=x(1) cmin=x(1) do 301 i=2,nxp ci=x(i) if(cmax.gt.ci) goto 302 cmax=ci 302 if(cmin.lt.ci) goto 301 cmin=ci 301 continue do 303 i=1,nxp 303 x(i)=(x(i)-cmin)/(cmax-cmin) endif ! ------------------------- Main Parameters Setting ------------------------- !- set Jianbo's parameter ie using variables that are controlled via ! calling arguments, and calculate rulers rl: ie = ie_value write(*,*) 'ie= ',ie if (ie.lt.0) goto 4001 rl(1)=1./2**ie goto 4002 4001 rl(1)=2**(-ie) 4002 write(*,*) ' rl(1)=', rl(1) do 113 i=2,9 113 rl(i)=rl(i-1)/2 !- set a few of Jianbo's parameters, using variables that are controlled ! via calling arguments: in = tan_motion_lim0 !- points for uncorrelation of tangential motion effects md = m !- embedding dimension mt = L !- delay time mstep = k(2)-k(1) !- step of k mmap = map !- number of k values if (md .le. 0) then write(*,*)'error--bad md value' stop endif !- calculate nn, the number of points for calculation. note that this is ! not the same as the number of points in a timeseries, because you don't ! necessarily need all those points to span the phase space. nn is set ! to the calculated value or 8000, whichever is smaller: nn = nxp - ( in + ((md-1) * mt) + (mstep * mmap) ) nn = MIN(nn, 8000) if (nn .le. 0) then write(*,*)'error--bad nn value' stop endif !- diagnostic messages: write(*,*) 'num pts in input series : nxp', nxp write(*,*) 'points for uncorrelation : in', in write(*,*) 'total points for calc. : nn', nn write(*,*) 'embedding dimension : md', md write(*,*) 'delay time : mt', mt write(*,*) 'step of k : mstep', mstep write(*,*) 'num of k values : mmap', mmap ! -------------------------- Main Computation Block ------------------------- do 3101 i=1,mmap 3101 mapped(i)=mstep*i do 3102 i=1,mmap do 3102 j=1,9 ssz(i,j)=0.0 nsz(i,j)=0 3102 continue nn1=nn+in do 1 i0=1,mt do 4 j0=i0+in,nn1 i=i0 j=j0 do 104 kk=1,md i1=i+(kk-1)*mt j1=j+(kk-1)*mt 104 a(map1,kk)=(x(i1)-x(j1))*(x(i1)-x(j1)) 339 d(map1)=0.0 do 105 kk=1,md 105 d(map1)=d(map1)+a(map1,kk) dd=d(map1) if(dd.le.0) goto 41 if(dd.ge.rl(1)) goto 41 il=1 if(dd.ge.rl(2)) goto 1000 il=2 if(dd.ge.rl(3)) goto 1000 il=3 if(dd.ge.rl(4)) goto 1000 il=4 if(dd.ge.rl(5)) goto 1000 il=5 if(dd.ge.rl(6)) goto 1000 il=6 if(dd.ge.rl(7)) goto 1000 il=7 if(dd.ge.rl(8)) goto 1000 il=8 if(dd.ge.rl(9)) goto 1000 il=9 1000 r(map1)=ALOG(dd) imp=1 9001 net=mapped(imp) do 1041 kk=1,md i1=i+(kk-1)*mt j1=j+(kk-1)*mt 1041 a(imp,kk)=(x(i1+net)-x(j1+net))*(x(i1+net)-x(j1+net)) d(imp)=0.0 do 1042 kk=1,md 1042 d(imp)=d(imp)+a(imp,kk) d2=d(imp) if(d2.le.0) goto 411 r(imp)=ALOG(d2) er=r(imp)-r(map1) ssz(imp,il)=ssz(imp,il)+er/2 nsz(imp,il)=nsz(imp,il)+1 411 imp=imp+1 if(imp.le.mmap) goto 9001 41 if(j.ge.nn1) goto 4 i=i+mt j=j+mt do 106 kk=1,md-1 106 a(map1,kk)=a(map1,kk+1) ldm=mt*(md-1) ii1=i+ldm jj1=j+ldm a(map1,md)=(x(ii1)-x(jj1))*(x(ii1)-x(jj1)) goto 339 4 continue 1 continue do 7001 il=1,9 do 7001 i=1,mmap mm=nsz(i,il) if(mm.eq.0) goto 7001 ssz(i,il)=ssz(i,il)/nsz(i,il) 7001 continue ! ------------------------ Prepare Output Variables ------------------------- do i=1,mmap out_log_disp(i,1) = ssz(i,1)-(ie+0)*ALOG(2.)/2. out_log_disp(i,2) = ssz(i,2)-(ie+1)*ALOG(2.)/2. out_log_disp(i,3) = ssz(i,3)-(ie+2)*ALOG(2.)/2. out_log_disp(i,4) = ssz(i,4)-(ie+3)*ALOG(2.)/2. out_log_disp(i,5) = ssz(i,5)-(ie+4)*ALOG(2.)/2. out_log_disp(i,6) = ssz(i,6)-(ie+5)*ALOG(2.)/2. out_log_disp(i,7) = ssz(i,7)-(ie+6)*ALOG(2.)/2. out_log_disp(i,8) = ssz(i,8)-(ie+7)*ALOG(2.)/2. out_log_disp(i,9) = ssz(i,9)-(ie+8)*ALOG(2.)/2. enddo out_lambda = ssz out_npts = nsz out_shells = SQRT(rl) end subroutine lambda_log_disp !======== end of file ========