2
\$\begingroup\$

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?

\$\endgroup\$
0
\$\begingroup\$

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)
```
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.