Improve `PowerSpectrum` estimator: add threhsold

[Datseris]

The power spectrum estimator will yield rather uncoverged results for small timeseries of regular signals, due to the noise induced by the FOurier transform affecting the signal a lot. We should add a threshold that reduces to 0 all spectral power less than this threshold. This would make the results much more reasonable for signals, and better for normalized. Here is the example code I have:

using DynamicalSystems
N1, N2, a = 101, 100001, 10

for N in (N1, N2)
    t = LinRange(0, 2*a*π, N)
    x = sin.(t) # periodic
    y = sin.(t .+ cos.(t/0.5)) # periodic, complex spectrum
    z = sin.(rand(1:15, N) ./ rand(1:10, N)) # random
    w = trajectory(Systems.lorenz(), N÷10; Δt = 0.1, Ttr = 100)[:, 1] # chaotic

    for q in (x, y, z, w)
        h = entropy(q, PowerSpectrum())
        n = entropy_normalized(q, PowerSpectrum())
        println("entropy: $(h), normalized: $(n).")
    end
end

[elavallee]

Hi, I would like to help with this issue. However, I am getting a different issue with the example code above. The line that is failing for me is this one:

w = trajectory(Systems.lorenz(), N÷10; Δt = 0.1, Ttr = 100)[:, 1] # chaotic

The error message is this:

MethodError: no method matching getindex(::Tuple{StateSpaceSets.StateSpaceSet{…}, StepRangeLen{…}}, ::Colon, ::Int64)
The function getindex exists, but no method is defined for this combination of argument types.

I am running Julia 1.12.2. I am new to contributing to open-source projects and to Julia, but if you point me in the right direction, I should be able to figure out the rest. It seems that Base.Tuple.getindex requires a new method that takes a StateSpaceSet as an argument. However, I am not sure how to fix that, and it may not be what is truly needed after all.

[kahaaga]

Hey @elavallee! Great that you want to contribute.

The error

The reason the example above gives you the error is that breaking changes has been introduced in subpackages of the DynamicalSystems.jl package since 2022.

The most recent documentation shows you how to fetch the relevant time series (hint: the time indices are returned by trajectory too, not just the actual system trajectory).

Check out the examples there - it should point you to how to generate the example time series.

The actual feature request

But the data input is really not that important here - we could use the PowerSpectrum estimator on any input time series.

What is needed is that the PowerSpectrum estimator gets a new keyword. This keyword will control whether or not some of the spectral power is set to 0 or not.

For some inspiration on how to do this, check of the https://github.com/JuliaDynamics/TimeseriesSurrogates.jl source code. We've used similar techniques there for some of the surrogate methods.

Feel free to give it a try, @elavallee, and don't hesitate to post any questions you may have here 👍

[elavallee]

Okay, I revised the code a bit to get it working. I see that entropy_normalized got deprecated, so I left it out. If I am not mistaken, Datseris' point is that the two evaluations of N1 and N2 should be closer numerically than they are. I will also look for examples in TimeseriesSurrogates.jl. Thanks @kahaaga.
`

using DynamicalSystems, ComplexityMeasures
N1, N2, a = 101, 100001, 10

for N in (N1, N2)
    println("N: $(N)")
    t = LinRange(0, 2*a*π, N)
    x = sin.(t) # periodic
    y = sin.(t .+ cos.(t/0.5)) # periodic, complex spectrum
    z = sin.(rand(1:15, N) ./ rand(1:10, N)) # random
    w, t = trajectory(Systems.lorenz(), N÷10; Δt = 0.1, Ttr = 100) # chaotic
    w    = w[:, 1]

    for q in (x, y, z, w)
        p, ~ = probabilities_and_outcomes(PowerSpectrum(), q)
        h = entropy(p)
        println("entropy: $(h).")
    end
end

[elavallee]

Hey @kahaaga,

I started a local change for this, and a question came up. I would like to know if the default value for the threshold should be 0.0 instead of 5.0. What is the more common case for users: having a low number of samples or a high number?
`

struct PowerSpectrum{T<:Real} <: OutcomeSpace
    δ::T

    function PowerSpectrum(δ::T = 5.0) where {T}
        new{T}(δ)
    end
end

[kahaaga]

Hey @kahaaga,

I started a local change for this, and a question came up. I would like to know if the default value for the threshold should be 0.0 instead of 5.0. What is the more common case for users: having a low number of samples or a high number?
`

struct PowerSpectrum{T<:Real} <: OutcomeSpace
    δ::T

    function PowerSpectrum(δ::T = 5.0) where {T}
        new{T}(δ)
    end
end

Good question!

The expected behavior according to the docstring is that no thresholding is performed (since thresholding is not mentioned). If we default to thresholding, it would be a breaking change. Therefore, the default behavior must be to not threshold. The docstring should clearly state this in the PR.

[kahaaga]

Hey @kahaaga,

I started a local change for this, and a question came up. I would like to know if the default value for the threshold should be 0.0 instead of 5.0. What is the more common case for users: having a low number of samples or a high number?
`

struct PowerSpectrum{T<:Real} <: OutcomeSpace
    δ::T

    function PowerSpectrum(δ::T = 5.0) where {T}
        new{T}(δ)
    end
end

In general, this is the philosophy of the package: the algorithms provided are completely agnostic about the input data. They should work - as specified in the original research papers - for all timeseries, regardless of how many points are in the timeseries.

We can provide the option to introduce parameters that change the behavior of the algorithm beyond standard behavior. For example: it can be helpful for some users to threshold if their timeseries are short. But this must be opt-in (with respect to the original description of the algorithm), not the default behaviour.

[elavallee]

Hey @kahaaga,

Thanks for the helpful explanation. It makes sense. I added a PR for this, but there is one thing I forgot to mention in the PR: the example of the random and chaotic systems does not lower their entropies with this change. Intuitively, that seems correct to me, but let me know what you think.