ANN: MaxMinFilters.jl - fast streaming maximum / minimum within moving window

Any references to those next-step refining stats?

Here is a minimal example of running mean, the same applies to higher order stats:

using Statistics

function runmean1(x::Vector{T}, n::Int64=10)::Vector{T} where {T<:Real}
    len = size(x,1)
    @assert n<len && n>1 "Argument n is out of bounds."
    out = zeros(len - n + 1)
    @inbounds for i = n:len
        out[i-n+1] = mean(view(x, i-n+1:i))
    end
    return out
end

function runmean2(x::Vector{T}, n::Int64=10)::Vector{T} where {T<:Real}
    len = size(x,1)
    @assert n<len && n>1 "Argument n is out of bounds."
    out = zeros(len - n + 1)
    sum::T = 0
    @inbounds for i = 1:n
        sum += x[i]
    end
    out[1] = sum / n
    @inbounds for i = n+1:size(x,1)
        sum += x[i] - x[i-n] # sum = previous sum - old point + new point
        out[i-n+1] = sum / n
    end
    return out
end

x = randn(1000000)
w = 100
@time m1 = runmean1(x, w)
@time m2 = runmean2(x, w)
m1 ≈ m2

For nan-stable version you can see my last PR to Indicators.jl:

Just FYI, you may be interested in transducers instead of stateful functions, e.g., if you want to squeeze out performance. See: https://github.com/tkf/Transducers.jl / https://clojure.org/reference/transducers

6 Likes

6 posts were split to a new topic: Using Transducers.jl instead of stateful functions for stream processing

Late to the party here. The algorithm in MaxMinFilters is the same one implemented in ImageFiltering; it just turned out the one in ImageFiltering had been optimized for small window sizes. A small tweak fixed the performance for large windows, so that now its performance is similar to MaxMinFilters.

ImageFiltering also supports delays, supports multiple dimensions, implements boundary conditions, and supports arbitrary functions. I think it’s a reasonable candidate for your “one-size-fits-all” package described in the OP, and I agree that we probably don’t need so many different packages that do the same thing.

6 Likes

Thanks, I have updated benchmark plot for ImageFiltering with that issue fixed - it has constant time now.

What delay are you talking about? I didn’t see any delay notations in ImageFiltering. Is it about filtering image sequence (video stream)?

Regarding to one general package, one impartant case is chunk-based stream processing for growing data with the ability to store filter state until the next data chunk arrives. I don’t think that applies to images, because you don’t usually have, say, 256x256 image that grows in its first dimension (512x256, 768x256, etc.).

1 Like

What delay are you talking about? I didn’t see any delay notations in ImageFiltering.

In the window specification you can supply a range, and that range doesn’t even have to include 0:

julia> mapwindow(extrema, 1:10, (-3:-1,), border="replicate")
10-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (1, 1)
 (1, 2)
 (1, 3)
 (2, 4)
 (3, 5)
 (4, 6)
 (5, 7)
 (6, 8)
 (7, 9)

one impartant case is chunk-based stream processing for growing data

True, but we routinely analyze images (movies) that are multiple terabytes where the last dimension is huge. Between Mmap and things like StreamingContainer I think we have almost everything one needs. I think the only missing component would be to split https://github.com/JuliaImages/ImageFiltering.jl/blob/5bcbe300a9af56b85a773b0edb156f72fa6f0594/src/mapwindow.jl#L408 so that it can take the Wedge W as input and return the W and cache upon output. That would allow you to preserve the state if you do need to make multiple calls to handle multiple chunks.

Maybe it’s already mentioned elsewhere, but I find it a bit of a discoverability issue that ImageFiltering contains many things that are useful outside image processing. For example, if I am looking for a function to process a vector, it’s tricky for me to think about ImageFiltering (unless I already know it).

1 Like

Yeah, we should probably rename it ArrayFiltering.jl or something. Would that fix the problem?

7 Likes

Yes, I think it’d make non-image processing people interested in the package more.

We’ve gotten to the point where an image is just an arbitrary-dimensional array where the eltype is <:ColorType. Other than that, JuliaImages is basically an extension of JuliaArrays, one where the notion of “index locality” really matters (i.e., entries with nearby indexes are in some sense related to one another). While some of the manipulations are pretty image-specific (e.g., ImageSegmentation, ImageFeatures), a lot of the material there can be applied to a much wider class of arrays. It makes it a little tricky to decide how to package things.

4 Likes

There is also some increased code complexity/readability, when you want to write generic functions for any dimensionality. And also, when you want to run along different dimension.

Sure, but we’ve worked quite hard to make sure that multidimensional code is easy to write, and very readable, in Julia (https://julialang.org/blog/2016/02/iteration). Is the small extra complexity of handling arbitrary indexes and arbitrary dimensionality really worth a whole bunch of packages that duplicate functionality?

2 Likes

Yeah, there’s a lot of overlap in DSP.jl’s conv / filt / resample functions - it would be good to do a pass through to see where they should be merged together.

1 Like

I agree. The onus should probably be on me to get that rolling: I’ve been talking about splitting/renaming ImageFiltering for some time, and the package is anyway due for an overhaul to take advantage of some of the nice features of modern Julia versions.

3 Likes

I think, when writing some new functions, it is easier to start with 1-dimensional case for simplicity, then debug it, and then it can be rewritten to multi-dimensional case, if applicable. And if there are some notable performance differences OR extra code complexity, we can leave both implementations.

It’s nice to have one implementation for 2 and more dimensions. But combining 1 and many has some extra burden, e.g. in DataStructures we have 1-dimensional CircularDeque, and there is no need to rewrite it for multidimensional deque with pushes along selected dims. If we really need that, it can be just added in separate source file.

Hi, recently there were thread started by @Kruxigt and I realised, that I need this functionality :slight_smile:

As I will need those moving max/min for a quantitative algorithms I am working now, I can spend some time on such consolidation or improving something in existing ones, cause I want to avoid creating one more package with this functionality.

This is what I am going to do, using transducers will allow operations on infinite streams.

@tim.holy, unfortunately my benchmarks are not reflecting that. Am I doing something wrong?

EDIT: now I see that it uses fast algorithm for extrema, but not minimum or maximum

julia> @btime ImageFiltering.mapwindow($maximum, $x, $1001);
  1.010 s (11021 allocations: 17.82 MiB)

julia> @btime movmax($x, $1001);
  9.996 ms (3 allocations: 7.64 MiB)

julia> length(x)
1000000
4 Likes

What functionality are you talking about?
There are different sets of possible functionality extensions:

  1. Steteless functions (to process small data chunks at once) => stateful functions for (realtime or larger than memory) data streams, that can be stopped and resumed later, so function state can be saved and loaded. This is also related to delays and boundary conditions handling.

  2. Special function types:

    2.1. Moving window functions => moving window with decimation, if window step > 1 (or “jumping” window)
    2.2. Special filters, like sorting filter for fast moving median and moving percentiles, etc.
    2.3. Regularly sampled data => irregularly sampled data, constant window length => variable window length.
    2.4. Event detectors: local maxima peak detectors, zero-crossing and threshold detectors, step detectors and so on. Also related to data formats for irregularly sampled events.

  3. Functions for 1-dimensional signals => multi-diomensional signals, like video streams.

  4. Composing chains of processing steps, like: “filter1 > decimate > sliding function > detect events > calc features > …” and so on. It is related to using transducers to feed data, or making functions and data buffers compatible for possible use with transducers. You may be interested in this thread:
    Using Transducers.jl instead of stateful functions for stream processing
    Speaking of transducers, I just don’t fully understand, how to use them for complex processing chains, if you need branching and joining data streams, select to store some intermediate data, log intermediate data for debug purpose, and so on.

All that extensions are partly compatible with each other, but not all.

Also, I think we should start with as simple use cases as possible.

2 Likes

Nice summary, I can see you were working a lot on this algorithm.

This is the most interesting for me - processing of very large or infinite streams from databases or other “larger than memory” sources. I don’t think stopping and resuming is crucial but probably you have more insight in this concept than me, so maybe just I cannot imagine a use case. There are coroutines, tasks and threading (see ThreadPools.jl), so they can handle interrupts in the data delivery quite nice.

Those functionalities sounds complicated, I am not sure the O(N) algorithm will handle them.

I guess we neeed an array of Deques to handle that. Probably it would be useful somewhere but I have only 1D cases to handle.

I don’t think they can handle that. Branching streams is not easy. If we have an iterator to database query stream it would be problematic to iterate it in 2 places with different cadence. The solution would be to open 2 independent connections and process them separately, eventually we can join those 2 processed streams later.

So for me simplest use case is to make this fast extrema algorithm work for any iterable things, not only AbstractArrays. I’ll try to prototype it soon, the discussion you linked is helpful. I did it long time ago in C# so this is feasible.

When I looked for an efficient moving/rolling/running median in Julia this was one of the first threads I found, but none of the mentioned packages could do a running median that was fast enough for my application.
So I took matters into my own hands and created one:

For window sizes above roughly 10 it is faster than sairus7’s https://github.com/sairus7/SortFilters.jl which can do arbitrary probability levels though. See the README.md for Benchmarks.

5 Likes