julia> x = rand(1000000);
julia> K = 9000;
julia> result2 = @btime movingmax($x, $K);
21.238 ms (4 allocations: 7.58 MiB)
julia> result3 = @btime rolling_max($x, $K);
24.017 ms (5 allocations: 7.64 MiB)
Yes they are close even for bigger N and K, I think the trick my algorithm is doing is no looking backward to the vector in each iteration, cause my Deque is keeping more information in a tuple. We have also a bit different memory access pattern (last/first in Deque). I have no idea why I did it this way. It is a rewrite from my old code in C#, but I recall I spent a lot of time on tuning this thing.
UPDATE: This package implements similar algorithm using custom ring buffers, it’s 2x faster than our approach.
Alright. With some help from @Jakub_Wronowski I fixed a bug in the code. I also switched to using a CircularDeque (this gives a big improvement in speed - more than 2x at least for some parameters). By doing that I could get close to the speed of the best libraries that are using a custom ring buffer. I will be satisfied with that for now!
Below is the current code:
function rolling_optima(vect::Any, win_size::Integer; debug=false)
roll_max = zeros(Float64, length(vect))
roll_min = zeros(Float64, length(vect))
opt_type = zeros(Int8, length(vect))
opt_index = zeros(Int64, length(vect))
d_max = CircularDeque{Int64}(win_size)
d_min = CircularDeque{Int64}(win_size)
win_forward = win_size
win_backward = win_forward*2
opt_ind = 0
#comp = 10
pushfirst!(d_max,1)
pushfirst!(d_min,1)
@inbounds for i in 2:length(vect)
if vect[i] > vect[i-1]
@inbounds for k in 1:length(d_max)
if vect[i] >= vect[first(d_max)]
popfirst!(d_max)
else break end
end
pushfirst!(d_max,i)
else
@inbounds for k in 1:length(d_min)
if vect[i] <= vect[first(d_min)]
popfirst!(d_min)
else break end
end
pushfirst!(d_min,i)
end
j = i + 1 - win_forward
if j >= 1
roll_max[i] = vect[last(d_max)]
roll_min[i] = vect[last(d_min)]
if j == last(d_max)
if roll_max[i] == roll_max[j]
if j-win_backward > 0
failed = false
for v in @view vect[j-win_backward:j-1]
if roll_max[j] < v
failed = true
break end
end
if !failed
opt_ind +=1
opt_index[opt_ind] = j
opt_type[opt_ind] = 1
end
end
end
pop!(d_max)
elseif j == last(d_min)
if roll_min[i] == roll_min[j]
if j-win_backward > 0
failed = false
for v in @view vect[j-win_backward:j-1]
if roll_min[j] > v
failed = true
break end
end
if !failed
opt_ind +=1
opt_index[opt_ind] = j
opt_type[opt_ind] = -1
end
end
end
pop!(d_min)
end
end
end
opt_index = opt_index[1:opt_ind]
opt_type = opt_type[1:opt_ind]
if debug
return opt_index, opt_type, roll_max, roll_min
else
return opt_index, opt_type
end
end
I attach the latest version of my code below. For many cases it almost matches the movmaxmin() from MaxMinFilters.jl which was the fastest implementation I could find. My code is surely not as good but adds the functionality of finding optima with the use of rolling maximum/minimum.
function rolling_optima(vect::Any, win_size::Integer; debug=false)
roll_max = Vector{Float64}(undef, length(vect))
roll_min = Vector{Float64}(undef, length(vect))
opt_type = Vector{Int8}(undef, length(vect))
opt_index = Vector{Int64}(undef, length(vect))
d_max = CircularDeque{Int64}(win_size)
d_min = CircularDeque{Int64}(win_size)
win_forward = win_size
win_backward = win_forward*2
opt_ind = 0
pushfirst!(d_max,1)
pushfirst!(d_min,1)
@inbounds for i in 2:length(vect)
if vect[i] > vect[i-1]
@inbounds for k in 1:length(d_max)
if vect[i] >= vect[first(d_max)]
popfirst!(d_max)
else break end
end
pushfirst!(d_max,i)
else
@inbounds for k in 1:length(d_min)
if vect[i] <= vect[first(d_min)]
popfirst!(d_min)
else break end
end
pushfirst!(d_min,i)
end
j = i + 1 - win_forward
if j >= 1
roll_max[i] = vect[last(d_max)]
roll_min[i] = vect[last(d_min)]
if j == last(d_max)
if roll_max[i] == roll_max[j]
if j-win_backward > 0
failed = false
# @inbounds for v in @view vect[j-win_backward:j-1]
check_back = @view vect[j-1:-1:j-win_backward]
@inbounds for k in 1:(win_backward)
if roll_max[j] < check_back[k]
failed = true
break end
end
if !failed
opt_ind +=1
opt_index[opt_ind] = j
opt_type[opt_ind] = 1
end
end
end
pop!(d_max)
elseif j == last(d_min)
if roll_min[i] == roll_min[j]
if j-win_backward > 0
failed = false
# @inbounds for v in @view vect[j-win_backward:j-1]
check_back = @view vect[j-1:-1:j-win_backward]
@inbounds for k in 1:(win_backward)
if roll_min[j] > check_back[k]
failed = true
break end
end
if !failed
opt_ind +=1
opt_index[opt_ind] = j
opt_type[opt_ind] = -1
end
end
end
pop!(d_min)
end
else
roll_max[i] = 0
roll_min[i] = 0
end
end
opt_index = opt_index[1:opt_ind]
opt_type = opt_type[1:opt_ind]
if debug
return opt_index, opt_type, roll_max, roll_min
else
return opt_index, opt_type
end
end
