Additional partial-randomization algorithms (#149)

* Update partial randomization to have a relative method

* Add spectrum-based randomization

* Fix surrogate plot but

* Add tests for partial randomization

* Tweak legend

* Tweaks to docstrings

* Switch to types for relative randomization

* Update docs for relative partial randomization

* Update docs

* Fix typos in docs for partial randomization

* Fix more typos

* Temporary commit

* Temporary commit

* Update plotting merge for 1.9

* Add partial randomization comparison figure

* Fix error

* Update for consistency

* More consistency

* Update tests for ArgumentErrors

* Update docs figure

* Make partial randomization defaults consistent

* Make partial randomization defaults consistent

* Update based on review

* Update docs with lorenz simulation

* Fix docs

* Update docstring for surrocompare

* Update for consistent style

* Fix unhidden doc example due to extra whitespace

* Update changelog and version number

* Fix version number in changelog to match most recent entry

* Consistent changelog versions

* Remove surrocompare

* No need to standardize for plots

---------

Co-authored-by: George Datseris <datseris.george@gmail.com>

Author: brendanjohnharris

.gitignore

@@ -7,4 +7,5 @@ Manifest.toml
 *.scss
 *.css
 .vscode
-*style.jl
+*style.jl
+*.code-workspace

CHANGELOG.md

@@ -1,7 +1,9 @@
 *Changelog is kept with respect to version 1.0. This software follows SymVer2.0*
 
-# 2.5.0
+# 2.6
+- Added surrogate methods: `RelativePartialRandomization`, `SpectralPartialRandomization`, `RelativePartialRandomizationAAFT`, and `SpectralPartialRandomizationAAFT`.
 
+# 2.5
 - Moved to Julia extensions (requiring julia v1.9).
 
 # 2.4
@@ -11,7 +13,8 @@
 - `pvalue` is now correctly overloaded from StatsAPI.jl.
 
 # 2.2
-- Implemented API for automating surrogate hypothesis tests using the new exported names `SurrogateTest` and `pvalue`.
+- Implemented API for automating surrogate hypothesis tests using the new exported names
+  `SurrogateTest` and `pvalue`.
 - New documentation section with an educative example of surrogate testing.
 
 # 2.1
@@ -29,17 +32,22 @@
 ## New features
 - New surrogate methods: `PartialRandomization`, `PartialRandomizationAAFT`, `TFTDAAFT`,
     `TFTDIAAFT`, and `RandomCascade`.
-- Using the `f` keyword, it is now possible to select whether circular shifting at
-    each wavelet coefficient level should be performed for `WLS` surrogates.
+- Using the `f` keyword, it is now possible to select whether circular shifting at each
+    wavelet coefficient level should be performed for `WLS` surrogates.
 
 # 1.3
 
 ## New features
-- The API now allows random number generator to be specified for all methods, enabling reproducibility
-- New surrogate methods: `PseudoPeriodicTwin`, `IrregularLombScargle` and `TFTDRandomFourier`.
+- The API now allows random number generator to be specified for all methods, enabling
+  reproducibility
+- New surrogate methods: `PseudoPeriodicTwin`, `IrregularLombScargle` and
+  `TFTDRandomFourier`.
 
 ## Bug fixes
-- Fixed error in `RandomFourier(true)` and `AAFT` surrogates, where phases were not correctly computed, leading to surrogates whose power spectra didn't match the power spectrum of the original signal as well as they should. No other Fourier-based surrogate were affected by this bug.
+- Fixed error in `RandomFourier(true)` and `AAFT` surrogates, where phases were not
+  correctly computed, leading to surrogates whose power spectra didn't match the power
+  spectrum of the original signal as well as they should. No other Fourier-based surrogate
+  were affected by this bug.
 
 # 1.2
 - New surrogate methods: `AutoRegressive`

Project.toml

@@ -2,7 +2,7 @@ name = "TimeseriesSurrogates"
 uuid = "c804724b-8c18-5caa-8579-6025a0767c70"
 authors = ["Kristian Agasøster Haaga <kahaaga@gmail.com>", "George Datseris"]
 repo = "https://github.com/JuliaDynamics/TimeseriesSurrogates.jl.git"
-version = "2.5.1"
+version = "2.6.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

docs/src/index.md

@@ -62,6 +62,10 @@ RandomFourier
 TFTDRandomFourier
 PartialRandomization
 PartialRandomizationAAFT
+RelativePartialRandomization
+RelativePartialRandomizationAAFT
+SpectralPartialRandomization
+SpectralPartialRandomizationAAFT
 AAFT
 TAAFT
 IAAFT
@@ -141,4 +145,4 @@ BiBTeX:
     title = {TimeseriesSurrogates.jl: a Julia package for generating surrogate data},
     journal = {Journal of Open Source Software}
 }
-```
+```

docs/src/methods/fourier_surrogates.md

@@ -26,12 +26,17 @@ surroplot(ts, s)
 ```
 
 
- ## Partial randomization
+## Partial randomization
 
- ### Without rescaling
+### Without rescaling
 
- [`PartialRandomization`](@ref) surrogates are similar to random phase surrogates,
- but allows for tuning the "degree" of phase randomization.
+[`PartialRandomization`](@ref) surrogates are similar to random phase surrogates, but allow for tuning the "degree" of phase randomization. 
+[`PartialRandomization`](@ref) use an algorithm introduced by [^Ortega1998], which draws random phases as:
+
+$$\phi \to \alpha \xi , \quad \xi \sim \mathcal{U}(0, 2\pi),$$
+
+where $\phi$ is a Fourier phase and $\mathcal{U}(0, 2\pi)$ is a uniform distribution.
+Tuning the randomization parameter, $\alpha$, produces a set of time series with varying degrees of randomness in their Fourier phases. 
 
 ```@example MAIN
 using TimeseriesSurrogates, CairoMakie
@@ -43,10 +48,65 @@ s = surrogate(ts, PartialRandomization(0.5))
 surroplot(ts, s)
 ```
 
+In addition to [`PartialRandomization`](@ref), we provide two other algorithms for producing partially randomized surrogates, outlined below.
+
+### Relative partial randomization
+
+The [`PartialRandomization`](@ref) algorithm corresponds to assigning entirely new phases to the Fourier spectrum with some degree of randomness, regardless of any deterministic structure in the original phases. As such, even for $\alpha = 0$ the surrogate time series can differ drastically from the original time series.
+
+By contrast, the [`RelativePartialRandomization`](@ref) procedure draws phases as:
+
+$$\phi \to \phi + \alpha \xi, \quad \xi \sim \mathcal{U}(0, 2\pi).$$
+
+With this algorithm, phases are progressively corrupted by higher values of $\alpha$: surrogates are identical to the original time series for $\alpha = 0$, equivalent to random noise for $\alpha = 1$, and retain some of the structure of the original time series when $0 < \alpha < 1$. This procedure is particularly useful for controlling the degree of chaoticity and non-linearity in surrogates of chaotic systems.
+
+### Spectral partial randomization
+
+Both of the algorithms above randomize phases at all frequency components to the same degree.
+To assess the contribution of different frequency components to the structure of a time series, the [`SpectralPartialRandomization`](@ref) algorithm only randomizes phases above a frequency threshold.
+The threshold is chosen as the lowest frequency at which the power spectrum of the original time series drops below a fraction $1-\alpha$ of its maximum value (such that the power contained above the frequency threshold is a proportion $\alpha$ of the total power, excluding the zero frequency).
+
+See the figure below for a comparison of the three partial randomization algorithms:
+```@example MAIN
+using DynamicalSystemsBase # hide
+function surrocompare(x, surr_types, params; color = ("#7143E0", 0.9), N=1000, linewidth=3, # hide
+                    transient=0, kwargs...) # hide
+    fig = Makie.Figure(resolution = (1080, 480), fontsize=22, kwargs...) # hide
+    for (j, a) in enumerate(surr_types) # hide
+        for (i, p) in enumerate(params) # hide
+            ax = Makie.Axis(fig[i,j]) # hide
+            hidedecorations!(ax) # hide
+            ax.ylabelvisible = true # hide
+            lines!(ax, surrogate(x, a(p...))[transient+1:transient+N]; color, linewidth) # hide
+            j == 1 && (ax.ylabel = "α = $(p)"; ax.ylabelfont = :bold) # hide
+            i == 1 && (ax.title = string(a)) # hide
+        end # hide
+    end # hide
+    colgap!(fig.layout, 30) # hide
+    rowgap!(fig.layout, 30) # hide
+    return fig # hide
+end # hide
+@inbounds function lorenz_rule!(du, u, p, t) # hide
+    σ = p[1]; ρ = p[2]; β = p[3] # hide
+    du[1] = σ*(u[2]-u[1]) # hide
+    du[2] = u[1]*(ρ-u[3]) - u[2] # hide
+    du[3] = u[1]*u[2] - β*u[3] # hide
+    return nothing # hide
+end # hide
+u0 = [0, 10.0, 0] # hide
+p0 = [10, 28, 8/3] # hide
+diffeq = (; abstol = 1e-9, reltol = 1e-9) # hide
+lorenz = CoupledODEs(lorenz_rule!, u0, p0; diffeq) # hide
+x = trajectory(lorenz, 1000; Ttr=500, Δt=0.025)[1][:, 1] # hide
+surr_types = [PartialRandomization, RelativePartialRandomization, SpectralPartialRandomization] # hide
+params = [0.0, 0.1, 0.25] # hide
+surrocompare(x, surr_types, params; transient=1000) # hide
+```
+
+
 ### With rescaling
 
-[`PartialRandomizationAAFT`](@ref) adds a rescaling step to the [`PartialRandomization`](@ref) surrogates to obtain surrogates that contain the same values as the original time
-series.
+[`PartialRandomizationAAFT`](@ref) adds a rescaling step to the [`PartialRandomization`](@ref) surrogates to obtain surrogates that contain the same values as the original time series. AAFT versions of [`RelativePartialRandomization`](@ref) and [`SpectralPartialRandomization`](@ref) are also available.
 
 ```@example MAIN
 using TimeseriesSurrogates, CairoMakie
@@ -57,7 +117,6 @@ s = surrogate(ts, PartialRandomizationAAFT(0.7))
 
 surroplot(ts, s)
 ```
-
 ## Amplitude adjusted Fourier transform (AAFT)
 
 

ext/TimeseriesSurrogatesUncertainData.jl

@@ -7,4 +7,4 @@ TimeseriesSurrogates.surrogate(x::UncertainDataset, method::Surrogate) = surroga
 TimeseriesSurrogates.surrogate(x::UncertainValueDataset, method::Surrogate) = surrogate(resample(x), method)
 TimeseriesSurrogates.surrogate(x::UncertainIndexDataset, method::Surrogate) = surrogate(resample(x), method)
 
-end
+end

ext/TimeseriesSurrogatesVisualizations.jl

@@ -3,27 +3,27 @@ module TimeseriesSurrogatesVisualizations
 using TimeseriesSurrogates, Makie
 
 function TimeseriesSurrogates.surroplot!(fig, x, a;
-        cx = "#191E44", cs = ("#7143E0", 0.9), nbins = 50, kwargs...
-    )
+    cx="#191E44", cs=("#7143E0", 0.9), nbins=50, kwargs...
+)
 
     t = 1:length(x)
     # make surrogate timeseries
     s = a isa Surrogate ? surrogate(x, method) : a
 
     # Time series
-    ax1, _ = Makie.lines(fig[1,1], t, x; color = cx, linewidth = 2)
-    Makie.lines!(ax1, t, s; color = cs, linewidth = 2)
+    ax1, _ = Makie.lines(fig[1, 1], t, x; color=cx, linewidth=2)
+    Makie.lines!(ax1, t, s; color=cs, linewidth=2)
 
     # Binned multitaper periodograms
     p, psurr = TimeseriesSurrogates.DSP.mt_pgram(x), TimeseriesSurrogates.DSP.mt_pgram(s)
-    ax3 = Makie.Axis(fig[2,1]; yscale = log10)
-    Makie.lines!(ax3, p.freq, p.power; color = cx, linewidth = 3)
-    Makie.lines!(ax3, psurr.freq, psurr.power; color = cs, linewidth = 3)
+    ax3 = Makie.Axis(fig[2, 1]; yscale=log10)
+    Makie.lines!(ax3, p.freq, p.power; color=cx, linewidth=3)
+    Makie.lines!(ax3, psurr.freq, psurr.power; color=cs, linewidth=3)
 
     # Histograms
-    ax4 = Makie.Axis(fig[3,1])
-    Makie.hist!(ax4, x; label = "original", bins = nbins, color = cx)
-    Makie.hist!(ax4, s; label = "surrogate", bins = nbins, color = cs)
+    ax4 = Makie.Axis(fig[3, 1])
+    Makie.hist!(ax4, x; label="original", bins=nbins, color=cx)
+    Makie.hist!(ax4, s; label="surrogate", bins=nbins, color=cs)
     Makie.axislegend(ax4)
 
     ax1.xlabel = "time step"
@@ -36,9 +36,9 @@ function TimeseriesSurrogates.surroplot!(fig, x, a;
 end
 
 function TimeseriesSurrogates.surroplot(x, s;
-    cx = "#191E44", cs = ("#7143E0", 0.9), nbins = 50, kwargs...)
-    fig = Makie.Figure(resolution = (500,500), kwargs...)
+    cx="#191E44", cs=("#7143E0", 0.9), nbins=50, kwargs...)
+    fig = Makie.Figure(resolution=(500, 500), kwargs...)
     surroplot!(fig, x, s; cx, cs, nbins)
 end
 
-end
+end

src/TimeseriesSurrogates.jl

@@ -31,6 +31,8 @@ include("methods/aaft.jl")
 include("methods/iaaft.jl")
 include("methods/truncated_fourier.jl")
 include("methods/partial_randomization.jl")
+include("methods/relative_partial_randomization.jl")
+include("methods/spectral_partial_randomization.jl")
 include("methods/wavelet_based.jl")
 include("methods/pseudoperiodic.jl")
 include("methods/pseudoperiodic_twin.jl")

src/methods/partial_randomization.jl

@@ -2,19 +2,21 @@ export PartialRandomization, PartialRandomizationAAFT
 
 """
     PartialRandomization(α = 0.5)
-`PartialRandomization` surrogates[^Ortega1998] are similar to [`RandomFourier`](@ref) phase 
-surrogates, but during the phase randomization step, instead of drawing phases from `[0, 2π]`,
-phases are drawn from `[0, 2π]*α`, where `α ∈ [0, 1]`. The authors refers to `α` as the 
-"degree" of phase randomization, where `α = 0` means `0 %` randomization and 
-`α = 1` means `100 %` randomization.
+`PartialRandomization` surrogates[^Ortega1998] are similar to [`RandomFourier`](@ref)
+phase surrogates, but during the phase randomization step, instead of drawing phases
+from `[0, 2π]`, phases are drawn from `[0, 2π]*α`, where `α ∈ [0, 1]`. The authors refers to `α` as the "degree" of phase randomization, where `α = 0` means `0 %` randomization
+and `α = 1` means `100 %` randomization.
+
+See [`RelativePartialRandomization`](@ref) and [`SpectralPartialRandomization`](@ref) for
+alternative partial-randomization algorithms
 
 [^Ortega1998]: Ortega, Guillermo J.; Louis, Enrique (1998). Smoothness Implies Determinism in Time Series: A Measure Based Approach. Physical Review Letters, 81(20), 4345–4348. doi:10.1103/PhysRevLett.81.4345
 """
 struct PartialRandomization{T} <: Surrogate
     α::T
 
-    function PartialRandomization(α::T) where T <: Real 
-        @assert 0 <= α <= 1
+    function PartialRandomization(α::T=0.5) where T <: Real
+        0 <= α <= 1 || throw(ArgumentError("α must be between 0 and 1"))
         return new{T}(α)
     end
 end
@@ -30,15 +32,15 @@ function surrogenerator(x::AbstractVector, rf::PartialRandomization, rng = Rando
     r = abs.(𝓕)
     ϕ = angle.(𝓕)
     coeffs = zero(r)
-    
-    init = (inverse = inverse, m = m, coeffs = coeffs, n = n, r = r, 
+
+    init = (inverse = inverse, m = m, coeffs = coeffs, n = n, r = r,
             ϕ = ϕ, shuffled𝓕 = shuffled𝓕)
     return SurrogateGenerator(rf, x, s, init, rng)
 end
 
 function (sg::SurrogateGenerator{<:PartialRandomization})()
-    inverse, m, coeffs, n, r, ϕ, shuffled𝓕 = 
-        getfield.(Ref(sg.init), 
+    inverse, m, coeffs, n, r, ϕ, shuffled𝓕 =
+        getfield.(Ref(sg.init),
         (:inverse, :m, :coeffs, :n, :r, :ϕ, :shuffled𝓕))
     s, rng = sg.s, sg.rng
     α = sg.method.α
@@ -52,20 +54,20 @@ end
 """
     PartialRandomizationAAFT(α = 0.5)
 
-`PartialRandomizationAAFF` surrogates are similar to [`PartialRandomization`](@ref) 
-surrogates[^Ortega1998], but adds a rescaling step, so that the surrogate has 
+`PartialRandomizationAAFF` surrogates are similar to [`PartialRandomization`](@ref)
+surrogates[^Ortega1998], but adds a rescaling step, so that the surrogate has
 the same values as the original time series (analogous to the rescaling done for
 [`AAFT`](@ref) surrogates).
-Partial randomization surrogates have, to the package authors' knowledge, not been 
+Partial randomization surrogates have, to the package authors' knowledge, not been
 published in scientific literature.
 
 [^Ortega1998]: Ortega, Guillermo J.; Louis, Enrique (1998). Smoothness Implies Determinism in Time Series: A Measure Based Approach. Physical Review Letters, 81(20), 4345–4348. doi:10.1103/PhysRevLett.81.4345
 """
 struct PartialRandomizationAAFT{T} <: Surrogate
     α::T
 
-    function PartialRandomizationAAFT(α::T) where T <: Real 
-        @assert 0 <= α <= 1
+    function PartialRandomizationAAFT(α::T=0.5) where T <: Real
+        0 <= α <= 1 || throw(ArgumentError("α must be between 0 and 1"))
         return new{T}(α)
     end
 end
@@ -92,4 +94,4 @@ function (sg::SurrogateGenerator{<:PartialRandomizationAAFT})()
     s[ix] .= x_sorted
 
     return s
-end
+end