Time-series search with early stopping

Question

I want to search for a pattern in a time-series while either ignoring the mean/shift/bias or the scale/standard deviation.

Consequently, I've written two functions.

The first function searches passes through the time-series, incrementally calculating the mean for each search-space sub-sequence and using this mean to normalize the sub-sequence before comparing it to a normalized query.

function euc_dist(data::Vector{Float64}, query::Vector{Float64}, current_best::Float64)::Float64
    sum = 0

    for (dd, qq) in zip(data, query)
        sum += (dd - qq) ^ 2
        if sum >= current_best
            break
        end
    end

    return sum
end


function run_ignore_bias(data::Vector{Float64}, query::Vector{Float64})::Tuple{Float64, Int}
    m = length(query)

    # normalize query in same manner data sub-sequence will be normalized
    query = query .- (sum(query) / m)

    current_best = Inf
    loc = -1

    # Keep current data in a double-size array to avoid using modulo
    # Basically, the data is stored twice and weird indexing arithmetic is used to avoid
    # using a LIFO queue and negative indexing.
    # Computational efficiency benefit unclear.
    t = zeros(Float64, 2*m)
    tz = zeros(Float64, m)

    run_sum = 0.

    for (d_i, dat) in enumerate(data)
        run_sum += dat

        t_idx = ((d_i - 1) % m) + 1
        t[t_idx] = dat
        t[t_idx + m] = dat

        if d_i >= m
            run_mean = run_sum / m

            # offset for search-space data
            s_off = (d_i % m) + 1
            # offset for search-space bound data
            s_bound_off = (d_i - 1) - (m - 1) + 1

            tz = t[s_off:s_off + m - 1] .- run_mean
            dist = euc_dist(tz, query, current_best)

            if dist < current_best
                current_best = dist
                loc = s_bound_off
            end

            run_sum -= t[s_off]
        end
    end

    return sqrt(current_best), loc
end

The second function does the same, except it normalizes according to the standard deviation.

function run_ignore_scale(data::Vector{Float64}, query::Vector{Float64})::Tuple{Float64, Int}
    m = length(query)

    # normalize scale query
    q_mean = sum(query) / m
    query = query / sqrt(sum(query.^2)/m - q_mean^2)
    current_best = Inf
    loc = -1

    # Keep current data in a double-size array to avoid using modulo
    # Basically, the data is stored twice and weird indexing arithmetic is used to avoid
    # using a LIFO queue and negative indexing.
    # Computational efficiency benefit unclear.
    t = zeros(Float64, 2*m)
    tz = zeros(Float64, m)

    run_sum = 0.
    run_sum2 = 0.

    for (d_i, dat) in enumerate(data)
        run_sum += dat
        run_sum2 += dat ^ 2

        t_idx = ((d_i - 1) % m) + 1
        t[t_idx] = dat
        t[t_idx + m] = dat

        if d_i >= m
            run_mean = run_sum / m
            # occasionally, a floating point error can cause this value to be negative, thus take the absolute value before sqrt
            run_std = sqrt(abs((run_sum2 / m) - (run_mean^2)))

            # offset for search-space data
            s_off = (d_i % m) + 1
            # offset for search-space bound data
            s_bound_off = (d_i - 1) - (m - 1) + 1

            tz = t[s_off:s_off + m - 1] / run_std
            dist = euc_dist(tz, query, current_best)
            @assert dist > 0

            if dist < current_best
                current_best = dist
                loc = s_bound_off
            end

            run_sum -= t[s_off]
            run_sum2 -= t[s_off] ^ 2
        end
    end

    return sqrt(current_best), loc
end

Here are the tests for both functions.

using Test

@testset "ignore bias" begin
    sig = [.2, .3, .5, -.4, .2, .3]
    data = vcat(zeros(2), sig .+ 1., zeros(8), 2*sig, zeros(4))
    val, idx = run_ignore_bias(data, sig)
    # should find shifted signal, but not scaled signal
    @test idx == 3
    @test isapprox(val, 0., atol=0.001)
end


@testset "ignore scale" begin
    sig = [.2, .3, .5, -.4, .2, .3]
    data = vcat(zeros(2), sig .+ 1., zeros(8), 2*sig, zeros(4))
    val, idx = run_ignore_scale(data, sig)
    # should find scaled signal, but not shifted
    @test idx == 17
    @test isapprox(val, 0., atol=0.001)
end


@testset "dist calc" begin
    dist = euc_dist([1., 2., 3.], [4., 5., 6.], Inf)
    @test isapprox(dist, 27.0, atol=0.001)

    dist = euc_dist([1., 2., 3.], [4., 5., 6.], 8.)
    @test isapprox(dist, 9.0, atol=0.001)
end

How do I reduce the code duplication between these two functions?

Answer 1

I created an iterator in Julia to remove the duplicated work of array viewing.

"""
Duplicate data and indexing arithmetic to avoid
using a LIFO queue or negative indexing.
"""
mutable struct t_iter
    data::Vector{Float64}
    t::Vector{Float64}
    length::Int
    q_len::Int

    function t_iter(data::Vector{Float64}, q_len::Int)
        return new(data, zeros(Float64, 2*length(data)), length(data), q_len)
    end
end


function Base.iterate(data::t_iter, d_i=1)
    if d_i >= data.length
        return nothing
    end

    dat = data.data[d_i]

    t_idx = ((d_i - 1) % data.q_len) + 1
    data.t[t_idx] = dat
    data.t[t_idx + data.q_len] = dat

    return ((d_i, dat, data.t), d_i+1)
end

The new functions and tests are as follows:

function run_ignore_bias(data::Vector{Float64}, query::Vector{Float64})::Tuple{Float64, Int}
    m = length(query)

    # normalize query in same manner data sub-sequence will be normalized
    query = query .- (sum(query) / m)

    current_best = Inf
    loc = -1

    tz = zeros(Float64, m)

    run_sum = 0.

    for (d_i, dat, t) in t_iter(data, m)
        run_sum += dat

        if d_i >= m
            run_mean = run_sum / m

            # offset for search-space data
            s_off = (d_i % m) + 1

            tz = t[s_off:s_off + m - 1] .- run_mean
            dist = euc_dist(tz, query, current_best)

            if dist < current_best
                current_best = dist
                loc = d_i - m + 1
            end

            run_sum -= t[s_off]
        end
    end

    return sqrt(current_best), loc
end


function run_ignore_scale(data::Vector{Float64}, query::Vector{Float64})::Tuple{Float64, Int}
    m = length(query)

    # normalize scale query
    q_mean = sum(query) / m
    query = query / sqrt(sum(query.^2)/m - q_mean^2)
    current_best = Inf
    loc = -1

    tz = zeros(Float64, m)

    run_sum = 0.
    run_sum2 = 0.

    for (d_i, dat, t) in t_iter(data, m)
        run_sum += dat
        run_sum2 += dat ^ 2

        if d_i >= m
            run_mean = run_sum / m
            # occasionally, a floating point error can cause this value to be negative, thus take the absolute value before sqrt
            run_std = sqrt(abs((run_sum2 / m) - (run_mean^2)))

            # offset for search-space data
            s_off = (d_i % m) + 1

            tz = t[s_off:s_off + m - 1] / run_std
            dist = euc_dist(tz, query, current_best)
            @assert dist > 0

            if dist < current_best
                current_best = dist
                loc = d_i - m + 1
            end

            run_sum -= t[s_off]
            run_sum2 -= t[s_off] ^ 2
        end
    end

    return sqrt(current_best), loc
end


sig = [.2, .3, .5, -.4, .2, .3]
data = vcat(zeros(2), sig .+ 1., zeros(8), 2*sig, zeros(4))
val, idx = run_ignore_bias(data, sig)
# should find shifted signal, but not scaled signal
@assert idx == 3
@assert isapprox(val, 0., atol=0.001)


sig = [.2, .3, .5, -.4, .2, .3]
data = vcat(zeros(2), sig .+ 1., zeros(8), 2*sig, zeros(4))
val, idx = run_ignore_scale(data, sig)
# should find scaled signal, but not shifted
@assert idx == 17
@assert isapprox(val, 0., atol=0.001)


dist = euc_dist([1., 2., 3.], [4., 5., 6.], Inf)
@assert isapprox(dist, 27.0, atol=0.001)

dist = euc_dist([1., 2., 3.], [4., 5., 6.], 8.)
@assert isapprox(dist, 9.0, atol=0.001)
```