Export fill_surrogate_test! (#161)

* Export fill_surrogate_test!

I hope the title is self-explanatory!

* put the inbounds elsewhere

* correct typo

* fix incorrect plot method

* replace resolution with size for figure

* fix columns typo

* hopefully fix the surrogate tests

* retrigger tests

* finally find and fix the doc bug

Author: Datseris

CHANGELOG.md

@@ -1,5 +1,9 @@
 *Changelog is kept with respect to version 1.0. This software follows SymVer2.0*
 
+# 2.7
+
+The function `fill_surrogate_test!` that was already in the source code and used by `pvalue` is now exported, so that if the users just want a parallelized estimation of surrogate discriminatory statistic they can just use this function.
+
 # 2.6
 - Added surrogate methods: `RelativePartialRandomization`, `SpectralPartialRandomization`, `RelativePartialRandomizationAAFT`, and `SpectralPartialRandomizationAAFT`.
 

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.6.4"
+version = "2.7.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

docs/src/collections/irregular_surrogates.md

@@ -10,12 +10,12 @@ This function uses the simulated annealing algorithm[^SchmitzSchreiber1999] to m
 ```@example MAIN
 using TimeseriesSurrogates, CairoMakie, Random
 
-# Example data: random AR1 process with a time axis with unevenly 
+# Example data: random AR1 process with a time axis with unevenly
 # spaced time steps
 rng = Random.MersenneTwister(1234)
 x = AR1(n_steps = 300)
 N = length(x)
-t = (1:N) - rand(N) 
+t = (1:N) - rand(N)
 
 # Use simulated annealing based on convergence of Lomb-Scargle periodograms
 # The time series is relatively long, so set tolerance a bit higher than default.

docs/src/index.md

@@ -33,6 +33,7 @@ surrogenerator
 
 ```@docs
 SurrogateTest
+fill_surrogate_test!
 pvalue(::SurrogateTest)
 ```
 
@@ -96,7 +97,7 @@ RandomCascade
 
 ### Other
 
-```@doc
+```@docs
 AutoRegressive
 ShuffleDimensions
 IrregularLombScargle

docs/src/methods/fourier_surrogates.md

@@ -30,13 +30,13 @@ surroplot(ts, s)
 
 ### Without rescaling
 
-[`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:
+[`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 Ortega et al., 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. 
+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

docs/src/methods/ppts.md

@@ -18,7 +18,7 @@ method = PseudoPeriodicTwin(d, τ, δ, ρ)
 surr_orbit = surrogate(x, method)
 
 # Get scalar surrogate time series from first and second column.
-s1, s2 = columns(surr_orbit)
+s1, s2 = surr_orbit[:, 1], surr_orbit[:, 2]
 
 # Scalar time series versus surrogate time series
 fig = Figure()

ext/TimeseriesSurrogatesVisualizations.jl

@@ -7,8 +7,8 @@ function TimeseriesSurrogates.surroplot!(fig, x, a;
 )
 
     t = 1:length(x)
-    # make surrogate timeseries
-    s = a isa Surrogate ? surrogate(x, method) : a
+    # make surrogate timeseries (allow possibility for `a` to be a timeseries or method)
+    s = a isa Surrogate ? surrogate(x, a) : a
 
     # Time series
     ax1, _ = Makie.lines(fig[1, 1], t, x; color=cx, linewidth=2)
@@ -37,7 +37,7 @@ end
 
 function TimeseriesSurrogates.surroplot(x, s;
     cx="#191E44", cs=("#7143E0", 0.9), nbins=50, kwargs...)
-    fig = Makie.Figure(resolution=(500, 500), kwargs...)
+    fig = Makie.Figure(size=(500, 500), kwargs...)
     surroplot!(fig, x, s; cx, cs, nbins)
 end
 

src/core/surrogate_test.jl

@@ -1,6 +1,6 @@
 using Random: AbstractRNG
 import StatsAPI: HypothesisTest, pvalue
-export SurrogateTest, pvalue
+export SurrogateTest, pvalue, fill_surrogate_test!
 using Base.Threads
 
 """
@@ -11,8 +11,11 @@ The tests requires as input a function `f` that given a timeseries (like `x`) it
 outputs a real number, and a method of how to generate surrogates.
 `f` is the function that computes the discriminatory statistic.
 
-Once called with [`pvalue`](@ref), the `test` stores the real value `rval` and surrogate
+Once called with [`pvalue`](@ref), the `test` estimates and then
+stores the real value `rval` and surrogate
 values `vals` of the discriminatory statistic in the fields `rval, vals` respectively.
+Alternatively, you can use [`fill_surrogate_test!`](@ref) directly if you don't care about
+the p-value.
 
 `SurrogateTest` automates the process described in the documentation page
 [Performing surrogate hypothesis tests](@ref).
@@ -24,7 +27,7 @@ values `vals` of the discriminatory statistic in the fields `rval, vals` respect
 - `rng = Random.default_rng()`: a random number generator.
 - `n::Int = 10_000`: how many surrogates to generate and compute `f` on.
 - `threaded = true`: Whether to parallelize looping over surrogate computations in
-  [`pvalue`](@ref) to the available threads (`Threads.nthreads()`).
+   to the available threads (`Threads.nthreads()`).
 """
 struct SurrogateTest{F<:Function, S<:SurrogateGenerator, X<:Real} <: HypothesisTest
     f::F
@@ -71,19 +74,28 @@ function Base.show(io::IO, ::MIME"text/plain", test::SurrogateTest)
     return
 end
 
-@inbounds function fill_surrogate_test!(test::SurrogateTest)
+"""
+    fill_surrogate_test!(test::SurrgateTest) → rval, vals
+
+Perform the computations foreseen by `test` and return
+the value of the discriminatory statistic for the real data `rval`
+and the distribution of values for the surrogates `vals`.
+
+This function is called by `pvalue`.
+"""
+function fill_surrogate_test!(test::SurrogateTest)
     if test.threaded
-        Threads.@threads for i in eachindex(test.vals)
+        @inbounds Threads.@threads for i in eachindex(test.vals)
             sgen = test.sgens[Threads.threadid()]
             test.vals[i] = test.f(sgen())
         end
     else
         sgen = first(sgens)
-        for i in eachindex(test.vals)
+        @inbounds for i in eachindex(test.vals)
             test.vals[i] = test.f(sgen())
         end
     end
-    return
+    return test.rval, test.vals
 end
 
 """
@@ -102,8 +114,7 @@ discriminatory statistic values. This is the case for statistics that quantify e
 For statistics that quantify autocorrelation, use `tail = :right` instead.
 """
 function pvalue(test::SurrogateTest; tail = :left)
-    fill_surrogate_test!(test)
-    (; rval, vals) = test
+    rval, vals = fill_surrogate_test!(test)
     if tail == :right
         p = count(v -> v ≥ rval, vals)
     elseif tail == :left

test/test_test.jl

@@ -13,6 +13,10 @@ test = SurrogateTest(q, x, RandomFourier(); n = 1000, rng)
 
 # the AR1 process is much more correlated than its surrogates!
 p = pvalue(test)
-p > 0.9
+@test p > 0.9
 p = pvalue(test; tail = :right)
-@test p < 0.1
+@test p < 0.1
+
+test = SurrogateTest(q, x, RandomFourier(); n = 1000, rng)
+rval, vals = fill_surrogate_test!(test)
+@test minimum(vals) ≤ rval ≤ maximum(vals)