Remove allocations from random fourier surrogates (#160)

* remove allocations from random fourier

* add correctness test for random fourier

* bump version

Author: Datseris

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

src/methods/randomfourier.jl

@@ -1,5 +1,9 @@
+using Random
+using LinearAlgebra
+
 export RandomFourier, FT
 
+
 """
     RandomFourier(phases = true)
 
@@ -36,26 +40,27 @@ function surrogenerator(x::AbstractVector, rf::RandomFourier, rng = Random.defau
     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{<:RandomFourier})()
-    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, phases = sg.s, sg.rng, sg.method.phases
 
     if phases
-        coeffs .= rand(rng, Uniform(0, 2π), n)
+        rand!(rng, Uniform(0, 2π), coeffs)
         shuffled𝓕 .= r .* exp.(coeffs .* 1im)
     else
         coeffs .= r .* rand(rng, Uniform(0, 2π), n)
         shuffled𝓕 .= coeffs .* exp.(ϕ .* 1im)
     end
-    s .= inverse * shuffled𝓕 .+ m
+    mul!(s, inverse, shuffled𝓕)
+    s .+= m
     return s
 end
 

test/all_method_tests.jl

@@ -3,7 +3,6 @@ using TimeseriesSurrogates
 using TimeseriesSurrogates.AbstractFFTs
 using TimeseriesSurrogates.Statistics
 using TimeseriesSurrogates.Random
-ENV["GKSwstype"] = "100"
 
 N = 500
 ts = cumsum(randn(N))
@@ -248,6 +247,12 @@ end
         s = surrogate(x, rf)
 
         @test length(s) == length(x)
+        # test that power spectrum is conserved
+        psx = abs2.(rfft(x))
+        pss = abs2.(rfft(s))
+        # For some reason I don't understand the last element of the spectrum
+        # is way off what is should be.
+        @test all(isapprox.(psx[1:end-1], pss[1:end-1]; atol = 1e-8))
     end
 
     @testset "random amplitudes" begin