From 158ab98f3ecb074bc13045f98bfb5e244ae13dad Mon Sep 17 00:00:00 2001 From: George Datseris Date: Mon, 17 Dec 2018 23:58:36 +0100 Subject: [PATCH] Remove ami/gmi code Closes #15 --- src/RecurrenceAnalysis.jl | 1 - src/delay.jl | 159 -------------------------------------- test/runtests.jl | 8 +- 3 files changed, 1 insertion(+), 167 deletions(-) delete mode 100644 src/delay.jl diff --git a/src/RecurrenceAnalysis.jl b/src/RecurrenceAnalysis.jl index b12280c..b91f81c 100644 --- a/src/RecurrenceAnalysis.jl +++ b/src/RecurrenceAnalysis.jl @@ -39,7 +39,6 @@ end include("matrices.jl") include("plot.jl") include("rqa.jl") -include("delay.jl") include("radius.jl") include("windowed.jl") diff --git a/src/delay.jl b/src/delay.jl deleted file mode 100644 index da75e02..0000000 --- a/src/delay.jl +++ /dev/null @@ -1,159 +0,0 @@ -""" - autocorrelation(x) - -Calculate autocorrelation of a vector for all positive lags. -""" -autocorrelation(x) = xcorr(x,x)[length(x):end] - -# Numbers of bin -# Freeman-Diaconis rule -function n_fd(x::AbstractVector) - first, last = extrema(x) - iqr = quantile(x,.75) - quantile(x,.25) - ceil(Integer, (last-first) / (2iqr * length(x)^(-1/3)) ) -end - -# Sturges -n_sturges(x::AbstractVector) = ceil(Integer, 1+log2(length(x))) - -function makebins(x, n) - x_ext = extrema(x) - range(x_ext[1], stop=x_ext[2], length=n+1) -end - -# Own simplified implementation of hist2d for square edges -function hist2!(h, x, edges) - n = length(edges)-1 - for r=1:n, c=1:n - h[r,c] = 0; - end - for i = 1:size(x,1) - r = searchsortedfirst(edges, x[i,1]) - 1 - c = searchsortedfirst(edges, x[i,2]) - 1 - if 1 <= r <= n && 1 <= c <= n - h[r,c] += 1 - end - end - h -end - -""" - ami(x, delay[, nbins]) - -Calculate the average mutual information of a signal for given delays. - -The delay can either be an integer (in which case the function returns an integer), -or a series of numbers (returning a vector) specified as an array of ascending integers, -a range, or a tuple with the limits of the range (minimum, maximum). - -The number of bins used to estimate the marginal and joint entropies of the -original and delayed signal can be given as a fixed integer, or as a -string that specifies a binning criterion. Currently available criteria are -Sturges (`"Sturges"`, the default value) or Freeman-Diaconis (`"FD"`). -""" -function ami(x, delay::Integer, nbins::Integer) - n = length(x)-delay - log2n = log2(n) - edges = makebins(x, nbins) - pxy = zeros(nbins,nbins) - hist2!(pxy, x[(1:n).+[0 delay]], edges) - px = sum(pxy, 2) - py = sum(pxy, 1) - ret = 0. - for ix=1:nbins, iy=1:nbins - if pxy[ix,iy] != 0 - pxpy = px[ix]*py[iy] - ret += pxy[ix,iy]/n*(log2(pxy[ix,iy]/pxpy)+log2n) - end - end - # Maximum entropy for a given nbins is log2(nbins) - ret/log2(nbins) -end - -function ami(x, delay::Union{Array, AbstractRange}, nbins::Integer) - d2 = maximum(delay) - miv = zeros(d2+1) - n = length(x)-d2 - log2n = log2(n) - # Pre-allocated arrays - pxy = zeros(nbins,nbins) - px = zeros(nbins,1) - py = zeros(1,nbins) - for d in delay - edges = makebins(x[1:n+d], nbins) - hist2!(pxy, x[(1:n).+[0 d]], edges) - sum!(px, pxy) - sum!(py, pxy) - for ix=1:nbins, iy=1:nbins - if pxy[ix,iy] != 0 - pxpy = px[ix]*py[iy] - miv[d+1] += pxy[ix,iy]/n*(log2(pxy[ix,iy]/pxpy)+log2n) - end - end - end - miv[1 .+ delay]./log2(nbins) -end - -ami(x, delay::Tuple{Integer,Integer}, nbins) = ami(x, delay[1]:delay[2], nbins) - -function ami(x, delay, nbins="Sturges") - # Define number of bins - dict_binfuns = Dict("Sturges" => n_sturges, - "FD" => n_fd) - !haskey(dict_binfuns,nbins) && error("incorrect definition of bins.") - binfun = dict_binfuns[nbins] - nbins = binfun(x) - ami(x, delay, nbins) -end - -""" - gmi(x, delay, radius) - -Calculate the generalized mutual information of a signal for given delays, based -on Rényi entropies. - -The delay can either be an integer (in which case the function returns an integer), -or a series of numbers (returning a vector) specified as an array of ascending integers, -a range, or a tuple with the limits of the range (minimum, maximum). - -The radius is a fixed value in the absolute scale of the signal, -used for the calculation of recurrence matrices to estimate -Rényi entropies (see `?recurrencematrix` for details). -Maximum norm and no scaling of distances are assumed. -""" -function gmi(x, delay::Integer, radius::Real) - rmfull = recurrencematrix(x,radius,scale=1) - n = length(x) - delay - # Entropy of the non-delayed series - rm1 = rmfull[1:n, 1:n] - h2 = -log2(nnz(rm1)/n^2) - # Joint entropy - rmj = rm1 .* rmfull[delay.+(1:n), delay.+(1:n)] - h2tau = -log2(nnz(rmj)/n^2) - # Maximum entropy for a given radius corresponds to uniform distribution - maxh2 = log2((maximum(x)-minimum(x))/radius) - (2h2 - h2tau)/maxh2 -end - -function gmi(x, delay::Union{Array, AbstractRange}, radius::Real) - d2 = maximum(delay) - rmfull = recurrencematrix(x,radius,scale=1) - h2tau = zeros(d2+1) - n = length(x) - d2 - # Entropy of the non-delayed series - rm1 = rmfull[1:n, 1:n] - h2 = -log2(nnz(rm1)/n^2) - # Joint entropy - rmj = spzeros(n, n) - for d in delay - rmj = rm1 .* rmfull[d.+(1:n), d.+(1:n)] - h2tau[d+1] = -log2(nnz(rmj)/n^2) - end - # Maximum entropy for a given radius corresponds to uniform distribution - maxh2 = log2((maximum(x)-minimum(x))/radius) - (2h2 .- h2tau[1 .+ delay])./maxh2 -end - -gmi(x, delay::Tuple{Integer,Integer}, radius) = gmi(x, delay[1]:delay[2], radius) - -# gmi(x::AbstractVector, delay) = gmi(x, delay, radius_mrr(x, .01)) diff --git a/test/runtests.jl b/test/runtests.jl index c7d7d4d..a393a2b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,13 +41,7 @@ sigma=10.; rho=28.; b=8/3; x0 = ones(3) lorenz_data = dynamical_system(x0, lorenz_eq(sigma,rho,b), .01, 1000) x = lorenz_data[501:2:end,1] -# Look for optimal delay -ami_def = ami(x, (1,12)) -@test findmin(ami_def)[2] == 9 -ami_fd = ami(x, (1,12), "FD") -@test findmin(ami_fd)[2] == 10 -ami_15 = ami(x, (1,12), 15) -gmi_10 = gmi(x, (1,12), 0.1) + # Look for optimal threshold dd, rr = sorteddistances(x, theiler=1) # Distance and recurrence matrices