New implementation for sifting

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/EmpiricalModeDecomposition.jl

@@ -9,7 +9,8 @@ using Base.Iterators
 
 import Base: iterate, IteratorSize, eltype
 
-export EMDIterable, SiftIterable
+export EMDIterable
+export sift, sift!
 export emd, eemd, ceemd, maketestdata
 
 include("testdata.jl")
@@ -33,7 +34,7 @@ function iterate(iter::EMDIterable, imf_prev=(iter.yvec,false))
     if imf_prev[2]
         return nothing
     elseif sum(abs, imf_prev[1]) > 0.0 && !ismonotonic(imf_prev[1])
-        imf = sift(imf_prev[1], iter.xvec, 0.1)
+        imf = sift(imf_prev[1], iter.xvec; reltol=0.1)
         return imf, (imf_prev[1]-imf, false)
     else
         return imf_prev[1], (imf_prev[1], true)
@@ -108,7 +109,7 @@ end
 Base.IteratorSize(::Type{CEEMDIterable{U,V,T}}) where {U,V,T} = Base.SizeUnknown()
 
 function Base.iterate(iter::CEEMDIterable)
-    imf0 = mean([sift(iter.yvec+noise,iter.xvec,0.1) for noise in iter.noise_ens])
+    imf0 = mean([sift(iter.yvec+noise, iter.xvec; reltol=0.1) for noise in iter.noise_ens])
     iter_ens = EMDIterable.(iter.noise_ens, [iter.xvec])
     state_ens = Tuple[]
     imf_state_ens = iterate.(iter_ens)
@@ -126,7 +127,7 @@ function Base.iterate(iter::CEEMDIterable, state::CEEMDState)
     elseif sum(abs,state.yvec)>vstop && !ismonotonic(state.yvec)
 
         # TODO: Format for clarity?
-        imf = vec(median(hcat([sift(state.yvec+noise[1], iter.xvec, 0.1) for noise in state.imf_state_ens]...),dims = 2))
+        imf = vec(median(hcat([sift(state.yvec+noise[1], iter.xvec; reltol=0.1) for noise in state.imf_state_ens]...),dims = 2))
 
         for iens in 1:length(state.iter_ens)
             r = iterate(state.iter_ens[iens], state.imf_state_ens[iens][2])
@@ -174,10 +175,10 @@ end
 
 function iaestimation(imf, xs)
     r = abs.(imf)
-    maxes = Int[]
+    maxs = Int[]
     mins = Int[]
-    localmaxmin!(imf, maxes, mins)
-    interpolate(xs[maxes], imf[maxes], xs, DierckXInterp())
+    localmaxmin!(imf, maxs, mins)
+    interpolate(xs[maxs], imf[maxs], xs, DierckXInterp())
 end
 
 

src/sift.jl

@@ -1,96 +1,135 @@
-"""
-Iterator for the sifting algorithm, with time series values `T` and positions `U`.
-"""
-struct SiftIterable{T<:AbstractVector, U<:AbstractVector}
-    "Values of the time series"
-    yvec::T
-    "Positions of the values"
-    xvec::U
-    "Number of steps of a stable sift after which the sifting is aborted"
-    stop_steps::Integer
-end
+using Printf
 
 """
-    SiftState
-
-Handle the intermediate results of the sifting.
+Iterator for the sifting algorithm.
 """
-mutable struct SiftState
-    "Values of the time series"
-    yvec
-    "Position of the time steps"
-    xvec
+mutable struct SiftIterable{solT <: AbstractVector, vecT <: AbstractVector, numT <: Real}
+    "Measurements/IMF"
+    y::solT
+    "Time series values"
+    x::vecT
     "Positions of the local maxima"
-    maxes::Vector{Int}
+    maxs::Vector{Int}
     "Positions of the local minima"
     mins::Vector{Int}
     "Indices of the zero crossings"
     crosses::Vector{Int}
-    "Sum of the abs value of yvec, used in the stopping criteria"
-    s
+    "Tolerance"
+    tol::numT
+    "Residuals used to check stopping criteria"
+    residual::numT
+    "Maximum number of iteration"
+    maxiter::Int
+    "Number of steps of a stable sift after which the sifting is aborted"
+    stop_steps::Int
     "Number of steps on which the number of zero crossings was fix"
-    fix_steps::Integer
+    fix_steps::Int
 end
 
-Base.IteratorSize(::Type{SiftIterable}) = Base.SizeUnknown()
+@inline converged(it::SiftIterable) = it.residual ≤ var(it.y) * it.tol
 
+@inline start(it::SiftIterable) = 0
 
-function iterate(iter::SiftIterable)
-    maxes = Int[]
-    mins = Int[]
-    s = sum(abs, iter.yvec)
-    @debug "Absolute sum of the data $s"
-    crosses = Int[]
+@inline done(it::SiftIterable, iteration::Int) = iteration ≥ it.maxiter || converged(it)
 
-    state = SiftState(iter.yvec, iter.xvec, maxes, mins, crosses, s, 0)
-    return state, state
-end
+# kernel sifting algorithm
+function iterate(it::SiftIterable, iteration::Int=start(it))
+    # Check for termination first
+    if done(it, iteration)
+        return nothing
+    end
 
+    localmaxmin!(it.y, it.maxs, it.mins)
+    maxlen = length(it.maxs)
+    minlen = length(it.mins)
 
+    it.fix_steps == it.stop_steps && return nothing
 
-function iterate(iter::SiftIterable, state::SiftState)
-    localmaxmin!(state.yvec, state.maxes, state.mins)
-    maxlen = length(state.maxes)
-    minlen = length(state.mins)
-    @debug "The number of minima: $minlen"
-    @debug "The number of maxima: $maxlen"
-    state.fix_steps == iter.stop_steps && return nothing
-    zerocrossing!(state.yvec, state.crosses)
-    abs(length(state.crosses) - maxlen - minlen) <= 1 && (state.fix_steps += 1)
+    zerocrossing!(it.y, it.crosses)
+
+    abs(length(it.crosses) - maxlen - minlen) <= 1 && (it.fix_steps += 1)
     if maxlen < 4 || minlen < 4
         return nothing
     end
-    smin = get_edgepoint(state.yvec, state.xvec, state.mins, first, isless)
-    smax = get_edgepoint(state.yvec, state.xvec, state.maxes, first, !isless)
-    emin = get_edgepoint(state.yvec, state.xvec, state.mins, last, isless)
-    emax = get_edgepoint(state.yvec, state.xvec, state.maxes, last, !isless)
+    smin = get_edgepoint(it.y, it.x, it.mins, first, isless)
+    smax = get_edgepoint(it.y, it.x, it.maxs, first, !isless)
+    emin = get_edgepoint(it.y, it.x, it.mins, last, isless)
+    emax = get_edgepoint(it.y, it.x, it.maxs, last, !isless)
+
     # TODO: Format for clarity?
-    maxTS = interpolate([first(state.xvec); state.xvec[state.maxes]; last(state.xvec)], [smax; state.yvec[state.maxes]; emax], state.xvec, DierckXInterp())
-    minTS = interpolate([first(state.xvec); state.xvec[state.mins]; last(state.xvec)], [smin; state.yvec[state.mins]; emin], state.xvec, DierckXInterp())
+    maxTS = interpolate([first(it.x); it.x[it.maxs]; last(it.x)], [smax; it.y[it.maxs]; emax], it.x, DierckXInterp())
+    minTS = interpolate([first(it.x); it.x[it.mins]; last(it.x)], [smin; it.y[it.mins]; emin], it.x, DierckXInterp())
     subs = 0.5*(maxTS + minTS)
-    state.s = sum(abs, subs)
-    @debug "Absolute sum of the not yet decomposed data $(state.s)"
-    state.yvec = state.yvec - subs
-    return state, state
+    it.residual = sum(abs, subs)
+
+    it.y .-= subs
+
+    # Return the residual at item and iteration number as state
+    it.residual, iteration + 1
 end
 
+# Utility functions
 
+function sift_iterator(y, x;
+                      reltol::Real = 0.1,
+                      maxiter::Int = 6,
+                      stop_steps::Int = 4)
+    maxs, mins = Int[], Int[]
+    crosses = Int[]
+    residual = sum(abs, y)
+    tolerance = 0.1
+    fix_steps = 0
+    # Return the iterable
+    return SiftIterable(y, x, maxs, mins, crosses,
+        tolerance, residual, maxiter, stop_steps, fix_steps)
+end
+
+"""
+    sift(y, x=1:length(y); kwargs...) -> imf
+
+Same as [`sift!`](@ref), but allocates a solution vector `imf` initialized with zeros.
 """
-    sift(y, xvec)
+sift(y, x=1:length(y); kwargs...) = sift!(zeros(eltype(y), size(y)), y, x; kwargs...)
 
-Sift the vector `y` whose points have x coordinates given by `xvec`.
 """
-function sift(yvec, xvec=1:length(yvec), tol=0.1)
-    ϵ = var(yvec) * tol
-    #@show ϵ
-    stop(state) = state.s <= ϵ
-    imf = nothing
-    num_steps = 0
-    for (i, step) in enumerate(halt(SiftIterable(yvec, xvec, 4), stop))
-        #@show sum(abs, step.yvec)
-        imf = step.yvec
-        num_steps = i
+    sift!(imf, y, x; kwargs...) -> imf
+
+# Arguments
+
+- `imf`: intrinsic mode function, will be updated in-place;
+- `y`: data measurements;
+- `x`: time stamp series.
+
+## Keywords
+
+- `reltol::Real = 0.1`: relative tolerance for the stopping condition;
+- `maxiter::Int = 10`: maximum number of iterations;
+- `verbose::Bool = false`: print method information;
+- `log::Bool = false`: keep track of the residual in each iteration.
+
+# Output
+
+**if `log` is `false`**
+
+- `imf`: intrinsic mode function.
+"""
+function sift!(imf, y, x;
+             reltol::Real = 0.1,
+             maxiter::Int = 6,
+             stop_steps::Int = 4,
+             log::Bool = false,
+             verbose::Bool = false,
+             kwargs...)
+
+    imf .= y
+    # Actually perform sifting
+    iterable = sift_iterator(imf, x; reltol, maxiter, stop_steps)
+
+    for (iteration, item) = enumerate(iterable)
+        verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual)
     end
-    #@show num_steps
-    imf
+
+    verbose && println()
+
+    iterable.y
 end

src/utils.jl

@@ -72,28 +72,4 @@ function interpolate(knotxvals::Vector, knotyvals::Vector, predictxvals::Abstrac
     spl = Dierckx.Spline1D(knotxvals, knotyvals, k=k)
     #@show spl, predictxvals
     Dierckx.evaluate(spl,predictxvals)
-end
-
-
-struct HaltingIterable{I, F}
-    iter::I
-    fun::F
-end
-
-function iterate(iter::HaltingIterable)
-    next = iterate(iter.iter)
-    return dispatch(iter, next)
-end
-
-function iterate(iter::HaltingIterable, (instruction, state))
-    if instruction == :halt return nothing end
-    next = iterate(iter.iter, state)
-    return dispatch(iter, next)
-end
-
-function dispatch(iter::HaltingIterable, next)
-    if next === nothing return nothing end
-    return next[1], (iter.fun(next[1]) ? :halt : :continue, next[2])
-end
-
-halt(iter::I, fun::F) where {I,F} = HaltingIterable{I,F}(iter, fun)
+end

test/runtests.jl

@@ -1,6 +1,5 @@
 using EmpiricalModeDecomposition
-import EmpiricalModeDecomposition: ismonotonic, localmaxmin!, get_edgepoint,
-    SiftIterable
+import EmpiricalModeDecomposition: ismonotonic, localmaxmin!, get_edgepoint
 using Test
 using Random
 
@@ -48,13 +47,10 @@ using Random
         @test get_edgepoint(y, t, maxes, last, !isless)  == 1.5
     end
 
-    @testset "SiftIterable" begin
+    @testset "sift" begin
         x = -1:0.1:2π+1
         measurements = sin.(x)
-        imf = zero(measurements)
-        for sift in Base.Iterators.take(SiftIterable(measurements,x,6),10)
-            imf = sift.yvec
-        end
+        imf = sift(measurements, x; maxiter=10, stop_steps=6)
         @test imf ≈ measurements
     end