Support for NaN input #21

wuhuangjianCN · 2023-03-22T01:39:07Z

There are many missing/corrupted data in real world cases, adding support for NaN in inputs would be great.

Firionus · 2023-03-23T11:14:54Z

EDIT: I accidentally sent this out too early, please excuse the edits to my draft.

Interesting suggestion. But I'm not sure what the expected outcome would be with these missing/NaN values.

Option 1: Poison Whole Window

Statistics just returns NaN/missing if just one of the values is NaN/missing:

using Statistics

julia> median([1,2,NaN])
NaN

julia> median([1,2,missing])
missing

Is this what you would expect? That a single NaN/missing value "poisons" the entire window? That is what RollingFunctions does (note they use asymmetric window tapering at the start and no tapering at the end):

julia> using RollingFunctions

julia> runmedian([1,2,3,4], 3)
4-element Vector{Float64}:
 1.0
 1.5
 2.0
 3.0

julia> runmedian([1,2,NaN,4], 3)
4-element Vector{Float64}:
   1.0
   1.5
 NaN
 NaN

julia> runmedian([missing,2,3,4], 3)
4-element Vector{Union{Missing, Float64}}:
  missing
  missing
  missing
 3.0

Option 2: Use Sorting Convention From Base

Base.sort puts NaN as if it was the highest number (see https://docs.julialang.org/en/v1/base/base/#Base.isless):

julia> sort([NaN,1,2,Inf,-Inf])
5-element Vector{Float64}:
 -Inf
   1.0
   2.0
  Inf
 NaN

I have not tested this exhaustively (there are MANY pitfalls with comparisons in the algorithm), but that might be what FastRunningMedian currently does:

julia> using FastRunningMedian

julia> running_median([1,2,NaN,4], 3, :asymmetric)
6-element Vector{Float64}:
   1.0
   1.5
   2.0
   4.0
 NaN
   4.0

However, there's no support for missing at the moment:

julia> running_median([missing,2,3,4], 3, :asymmetric)
ERROR: MethodError: no method matching running_median(::Vector{Union{Missing, Int64}}, ::Int64, ::Symbol)

Closest candidates are:
  running_median(::AbstractVector{T}, ::Integer, ::Any) where T<:Real
   @ FastRunningMedian ~/.julia/packages/FastRunningMedian/TyoOP/src/FastRunningMedian.jl:28
  running_median(::AbstractVector{T}, ::Integer) where T<:Real
   @ FastRunningMedian ~/.julia/packages/FastRunningMedian/TyoOP/src/FastRunningMedian.jl:28

Stacktrace:
 [1] top-level scope
   @ REPL[78]:1

Option 3: Use Inverse Sorting Convention From Base

This is what you get with

julia> sort([NaN,1,2,Inf,-Inf], lt=<)
5-element Vector{Float64}:
 NaN
 -Inf
   1.0
   2.0
  Inf

I think this might be used by SortFilters (again, not tested exhaustively):

julia> using SortFilters

julia> movsort([1,2,3,4], 3, .5)
4-element Vector{Float64}:
 1.0
 1.0
 2.0
 3.0

julia> movsort([1,2,NaN,4], 3, .5)
4-element Vector{Float64}:
 1.0
 1.0
 1.0
 2.0

julia> movsort([4,3,NaN,1], 3, .5)
4-element Vector{Float64}:
   4.0
   4.0
   3.0
 NaN

julia> movsort([missing,2,3,4], 3, .5)
ERROR: TypeError: non-boolean (Missing) used in boolean context
Stacktrace:
 [1] movsort!(out::Vector{Float64}, in::Vector{Union{Missing, Int64}}, window::Int64, p::Float64)
   @ SortFilters ~/.julia/packages/SortFilters/lSAIo/src/SortFilters.jl:39
 [2] movsort(in::Vector{Union{Missing, Int64}}, window::Int64, p::Float64)
   @ SortFilters ~/.julia/packages/SortFilters/lSAIo/src/SortFilters.jl:17
 [3] top-level scope
   @ REPL[56]:1

Option 4: Ignore NaN/missing Unless They are the Only Thing in Window

There's another approach, which would be to skip all NaN/missing values if there's another one in the window and return NaN/missing otherwise:

julia> nonexistent_running_median([1,2,3,4], 3, :asymmetric)
6-element Vector{Float64}:
 1.0
 1.5
 2.0
 3.0
 3.5
 4.0

julia> nonexistent_running_median([1,2,NaN,4], 3, :asymmetric)
6-element Vector{Float64}:
 1.0
 1.5
 1.5
 3.0
 4.0
 4.0

julia> nonexistent_running_median([missing,2,3,4], 3, :asymmetric)
6-element Vector{Union{Missing, Float64}}:
 missing
 2.0
 2.5
 3.0
 3.5
 4.0

The behavior here is that the window still uses the elements from its supporting points but the actual number of used elements varies depending on how many NaN/missing are in there.

In consequence, with this option I would also expect:

julia> nonexistent_running_median([missing,NaN,3,4], 3, :asymmetric)
6-element Vector{Float64}:
 missing
 missing
 3.0
 3.5
 3.5
 4.0

Conclusion

For me, Option 4 sounds the most promising right now, though also probably the hardest to implement.

Option 2/3 I dislike because whether to put NaN as lowest or highest number is a completely arbitrary choice.

Option 1 is the only one that makes any mathematical sense, but is probably less useful in practice than the other options, even if the other options introduce subtle errors (Option 2 biases the median upwards, Option 3 biases the median downwards and Option 3 biases the median left/right depending on whether the missing value is to the right/left of the window center).

Which option would you like?

wuhuangjianCN · 2023-03-24T03:22:04Z

Thanks for your response. Option 4 is what I need.
Looking forward to your updates~

wuhuangjianCN · 2023-03-27T01:11:54Z

EDIT: I accidentally sent this out too early, please excuse the edits to my draft.

Interesting suggestion. But I'm not sure what the expected outcome would be with these missing/NaN values.

Option 1: Poison Whole Window

Statistics just returns NaN/missing if just one of the values is NaN/missing:
using Statistics

julia> median([1,2,NaN])
NaN

julia> median([1,2,missing])
missing
Is this what you would expect? That a single NaN/missing value "poisons" the entire window? That is what RollingFunctions does (note they use asymmetric window tapering at the start and no tapering at the end):
julia> using RollingFunctions

julia> runmedian([1,2,3,4], 3)
4-element Vector{Float64}:
 1.0
 1.5
 2.0
 3.0

julia> runmedian([1,2,NaN,4], 3)
4-element Vector{Float64}:
   1.0
   1.5
 NaN
 NaN

julia> runmedian([missing,2,3,4], 3)
4-element Vector{Union{Missing, Float64}}:
  missing
  missing
  missing
 3.0
Option 2: Use Sorting Convention From Base

Base.sort puts NaN as if it was the highest number (see https://docs.julialang.org/en/v1/base/base/#Base.isless):
julia> sort([NaN,1,2,Inf,-Inf])
5-element Vector{Float64}:
 -Inf
   1.0
   2.0
  Inf
 NaN
I have not tested this exhaustively (there are MANY pitfalls with comparisons in the algorithm), but that might be what FastRunningMedian currently does:
julia> using FastRunningMedian

julia> running_median([1,2,NaN,4], 3, :asymmetric)
6-element Vector{Float64}:
   1.0
   1.5
   2.0
   4.0
 NaN
   4.0
However, there's no support for missing at the moment:
julia> running_median([missing,2,3,4], 3, :asymmetric)
ERROR: MethodError: no method matching running_median(::Vector{Union{Missing, Int64}}, ::Int64, ::Symbol)

Closest candidates are:
  running_median(::AbstractVector{T}, ::Integer, ::Any) where T<:Real
   @ FastRunningMedian ~/.julia/packages/FastRunningMedian/TyoOP/src/FastRunningMedian.jl:28
  running_median(::AbstractVector{T}, ::Integer) where T<:Real
   @ FastRunningMedian ~/.julia/packages/FastRunningMedian/TyoOP/src/FastRunningMedian.jl:28

Stacktrace:
 [1] top-level scope
   @ REPL[78]:1
Option 3: Use Inverse Sorting Convention From Base

This is what you get with
julia> sort([NaN,1,2,Inf,-Inf], lt=<)
5-element Vector{Float64}:
 NaN
 -Inf
   1.0
   2.0
  Inf
I think this might be used by SortFilters (again, not tested exhaustively):
julia> using SortFilters

julia> movsort([1,2,3,4], 3, .5)
4-element Vector{Float64}:
 1.0
 1.0
 2.0
 3.0

julia> movsort([1,2,NaN,4], 3, .5)
4-element Vector{Float64}:
 1.0
 1.0
 1.0
 2.0

julia> movsort([4,3,NaN,1], 3, .5)
4-element Vector{Float64}:
   4.0
   4.0
   3.0
 NaN

julia> movsort([missing,2,3,4], 3, .5)
ERROR: TypeError: non-boolean (Missing) used in boolean context
Stacktrace:
 [1] movsort!(out::Vector{Float64}, in::Vector{Union{Missing, Int64}}, window::Int64, p::Float64)
   @ SortFilters ~/.julia/packages/SortFilters/lSAIo/src/SortFilters.jl:39
 [2] movsort(in::Vector{Union{Missing, Int64}}, window::Int64, p::Float64)
   @ SortFilters ~/.julia/packages/SortFilters/lSAIo/src/SortFilters.jl:17
 [3] top-level scope
   @ REPL[56]:1
Option 4: Ignore NaN/missing Unless They are the Only Thing in Window

There's another approach, which would be to skip all NaN/missing values if there's another one in the window and return NaN/missing otherwise:
julia> nonexistent_running_median([1,2,3,4], 3, :asymmetric)
6-element Vector{Float64}:
 1.0
 1.5
 2.0
 3.0
 3.5
 4.0

julia> nonexistent_running_median([1,2,NaN,4], 3, :asymmetric)
6-element Vector{Float64}:
 1.0
 1.5
 1.5
 3.0
 4.0
 4.0

julia> nonexistent_running_median([missing,2,3,4], 3, :asymmetric)
6-element Vector{Union{Missing, Float64}}:
 missing
 2.0
 2.5
 3.0
 3.5
 4.0
The behavior here is that the window still uses the elements from its supporting points but the actual number of used elements varies depending on how many NaN/missing are in there.

In consequence, with this option I would also expect:
julia> nonexistent_running_median([missing,NaN,3,4], 3, :asymmetric)
6-element Vector{Float64}:
 missing
 missing
 3.0
 3.5
 3.5
 4.0
Conclusion

For me, Option 4 sounds the most promising right now, though also probably the hardest to implement.

Option 2/3 I dislike because whether to put NaN as lowest or highest number is a completely arbitrary choice.

Option 1 is the only one that makes any mathematical sense, but is probably less useful in practice than the other options, even if the other options introduce subtle errors (Option 2 biases the median upwards, Option 3 biases the median downwards and Option 3 biases the median left/right depending on whether the missing value is to the right/left of the window center).

Which option would you like?

I wrote a simple implementation for Option 4:


function medfilter(inData;edgeLen=24*2)
    # running_median that ignore nan
    q_nan=isnan.(inData)
    x=copy(inData)
    x[q_nan] .= -Inf
    a=running_median(x, edgeLen * 2 + 1)
    x[q_nan] .= Inf
    b=running_median(x, edgeLen * 2 + 1)
    # window without valid data will return NaN, because Inf - Inf = NaN
    return (a + b)./2
end

However, the performance can be improved and NaN will return when there are more than half of NaN inside a window.

a=5 .*sin.(LinRange(0,20*pi,1000))+rand(1000)
a[14:27] .= -Inf
a[:4] .= -Inf
b= running_median(a, 5 * 2 + 1)
using PyPlot
plt=get_plt()
plt.clf()
plt.plot(a,labels="input")
plt.plot(b,labels="output")

plt.xlim([0,50])
plt.gcf()

Firionus · 2023-03-29T19:09:30Z

Thanks for your response. Option 4 is what I need.
Looking forward to your updates~

In an earlier version of your answer (I got it via mail) you mentioned the behavior of MATLAB. I looked at it and happen to really like it: https://www.mathworks.com/help/matlab/ref/movmedian.html.
By default, they use Option 1, so if there is just 1 NaN in the window, the median is also NaN. This makes a lot of sense as default behavior since it is mathematically the most sound option.
However, MATLAB also allows you to set the flag "omitnan" and get the behavior of Option 4.

I wrote a simple implementation for Option 4:

That's elegant, but not quite right. Inserting Inf or -Inf pushes the median in those windows up or down by some amount that strongly depends on the distribution of the data, and is not guaranteed to average out afterwards. It is not the same as completely omitting the NaN.

Anyway, I'd like to copy MATLAB's behavior. Here's the steps I'd like to take:

Write a slow reference implementation for Option 1 and 4 with RollingFunctions (I'll post them here shortly).
Write fuzzing tests for Option 1 and 4 based on the reference implementation
Implement Option 1 (default)
Implement Option 4
Documentation
Benchmarks

The biggest hurdle will be an elegant implementation that works well with both options. Option 1 will probably want to track the index of the newest NaN in the window, and output NaN until it is out again. Option 4 will likely need the ability to roll with less elements than the maximum window size, which currently is impossible.

I'll start with the reference implementation and I'll also do the Fuzz tests in the near future.

I'll also split out an issue regarding more data types, like Union{Float64, Missing} or Dates.

Firionus · 2023-03-29T20:01:58Z

Alright, reference implementation for Option 4 is done:

julia> using RollingFunctions, Statistics

julia> function median_ignorenan(x)
         mask = @. !isnan(x)
         masked = @view x[mask]
         masked |> isempty && return NaN
         median(masked)
       end
median_ignorenan (generic function with 1 method)

julia> running(median_ignorenan, [NaN,2,3,4], 3)
4-element Vector{Float64}:
 NaN
   2.0
   2.5
   3.0

julia> running(median_ignorenan, [1,2,NaN,4], 3)
4-element Vector{Float64}:
 1.0
 1.5
 1.5
 3.0

Until FastRunningMedian supports it, the above (maybe replace running with rolling, depending on the boundary conditions you want - check the docs of RollingFunctions please!) should probably be your workaround for NaNs, even though it is less performant than this library. If you want Option 1, just use regular RollingFunctions.

wuhuangjianCN · 2023-03-30T02:48:27Z

That's elegant, but not quite right. Inserting Inf or -Inf pushes the median in those windows up or down by some amount that strongly depends on the distribution of the data, and is not guaranteed to average out afterwards. It is not the same as completely omitting the NaN.

Yes, that approah misbehaves when encounter too many NaNs inside a running window. Yet, I still use it for now because the performance of running median is critical in my task, and your package is faster than calling the pandas package by times.

Good luck with your future update~

berjine · 2023-06-19T13:17:05Z

Hi, could it work to just replace the median calls by nanmedian from NaNStatistics ?
Because the implementation using RollingFunctions seems much slower than the standard one from FastRunningMedian

Firionus · 2023-06-19T19:34:15Z

Hi @berjine 👋
you can benchmark whether ‘nanmedian’ combined with RollingFunctions is faster than what I provided as a workaround above. Feel free to share your results here. I wouldn’t expect big improvements though.

Unfortunately I haven’t had much time for programming the NaN handling. The testing is in place on the ‘master’ branch though, so if someone wants to give the implementation a go, feel free to do so and put in a PR. If you start working on it, let us know in this issue.

Last time I started working on this I remember being so fed up with modifying my complicated algorithm that I explored SkipLists some more instead 😂. I’ll see whether I can carve out some more time for this issue in the future, but it’s a hobby project, so no guarantees.

berjine · 2023-06-20T09:05:19Z

Sure !

A is running(median_ignorenan, test_data, windowsize)
B is running(nanmedian, test_data, windowsize)
C is running_median(test_data, windowsize), added to see up to where it's possible to go !

With 1M points, 1% random NaN and 1001-points window :
9.557 s (6999732 allocations: 26.79 GiB) for A
11.550 s (1000007 allocations: 7.64 GiB) for B
69.380 ms (29 allocations: 7.70 MiB) for C, 138x faster than A (!!!)

Same data but 101-points window :
2.075 s (5999986 allocations: 2.61 GiB) for A
3.655 s (1000006 allocations: 869.71 MiB) for B
59.894 ms (25 allocations: 7.65 MiB) for C, 35x faster than A

With 10k points, 1% random NaN and 101-points window :
12.643 ms (59986 allocations: 26.58 MiB) for A
33.365 ms (10006 allocations: 8.66 MiB) for B
547.900 μs (25 allocations: 94.91 KiB) for C, 23x faster than A

It confirms your workaround works great, and that FastRunningMedian is fantastic haha

A possible workaround is to interpolate the NaNs to suppress them beforehand
This step take 96ms for 1M points at 1% NaNs with SteffenMonotonicInterpolation from Interpolations.jl

This works quite well I think, as long as there are no NaNs at the beginning or end of the data vector x)

We can note that both A and B implementations seem to be shifted by half the window length, I don't know why though

Firionus · 2023-06-20T10:24:31Z

Thanks for the comparison. Interpolating NaNs might be the best workaround yet in terms of performance 👍

We can note that both A and B implementations seem to be shifted by half the window length, I don't know why though

That's because RollingFunctions.jl uses different boundary conditions than FastRunningMedian.jl. If I remeber correctly (untested code ahead), you can match the two with either:

running_median(in, N, tapering=:asymmetric)[1:length(in)] and running(median_ignorenan, in, N)
running_median(in, N, tapering=:none) and rolling(median_ignorenan, in, N)

Consult the Excel sheet in the README or the docstrings for a detailed overview of the boundary conditions in this library.

Firionus · 2023-06-20T14:07:08Z

I got around to the implementation of Option 1 (default) and Option 4 today. 🎉

It is available on the master branch. You can install it with ]add https://github.com/Firionus/FastRunningMedian.jl.git#master. Ignore the NaNs like this:

running_median(in, w, nan=:ignore)

Default behavior is nan=:include, where a single NaN will turn the entire window NaN.

I'd appreciate if both of you could give it a try and compare it to your existing solutions. Also, if you have Feedback on Performance or the API, please feel free to write what's on your mind.

berjine · 2023-06-20T15:01:20Z

Fantastic ! It works with the same awesome performance !
79.248 ms (5468 allocations: 8.09 MiB) with nan=:ignore on 1M points, 1% NaN and 1001-points window

Maybe you would want to detect the presence of NaN or missing values and autoselect the nan=:ignore mode ?

Firionus · 2023-06-20T18:19:49Z

Maybe you would want to detect the presence of NaN or missing values and autoselect the nan=:ignore mode ?

See above:

Option 1 (nan=:include) is the only one that makes any mathematical sense, but is probably less useful in practice than the other options, even if the other options introduce subtle errors (Option 2 biases the median upwards, Option 3 biases the median downwards and Option 4 (nan=:ignore) biases the median left/right depending on whether the missing value is to the right/left of the window center).

MATLAB also does it like this and I‘m pretty convinced it is the right approach. That‘s why nan=:include is the default.

Firionus · 2023-06-25T10:17:07Z

Very quick benchmark comparison of how much the NaN handling slowed down the performance:

Before (0.2.0):

julia> @benchmark running_median($x, 101)
BenchmarkTools.Trial: 88 samples with 1 evaluation.
 Range (min … max):  54.663 ms … 67.524 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     56.883 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   57.382 ms ±  1.701 ms  ┊ GC (mean ± σ):  0.19% ± 0.47%

            ▁▁▂█▄ ▁▁                                           
  ▃▁▁▁▁▃▃▃▁▅████████▅▅█▅▃▅█▁▃▅▃▅▁▆▁▁▅▃▃▃▁▁▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  54.7 ms         Histogram: frequency by time        62.8 ms <

 Memory estimate: 7.65 MiB, allocs estimate: 25.

After (9b17e62):

julia> @benchmark running_median($x, 101)
BenchmarkTools.Trial: 85 samples with 1 evaluation.
 Range (min … max):  57.673 ms … 72.718 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     58.444 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   58.995 ms ±  1.753 ms  ┊ GC (mean ± σ):  0.17% ± 0.43%

         ▂█   ▂                                                
  ▄█▇▄▁▇▅██▇▇▅█▄▄▁▁▄▁▁█▅▅▄▅▅▅▄▁▅▅▁▁▁▁▄▁▄▁▄▇▁▁▄▇▄▁▁▁▄▁▅▇▄▁▄▁▁▄ ▁
  57.7 ms         Histogram: frequency by time        60.9 ms <

 Memory estimate: 7.65 MiB, allocs estimate: 64.

That's a slowdown of 3 %. It's not great but NaNs need to be handled for correctness. I'll ship it now and we can include possible performance improvements in the next patch release.

You should soon be able to update FastRunningMedian to v0.2.1 and get the new NaN handling feature.

As of v0.2.1, FastRunningMedian doesn't output garbage with NaNs anymore 🎉 (tracking issue: Firionus/FastRunningMedian.jl#21). I've adjusted the README slightly to include this information.

Firionus mentioned this issue Mar 29, 2023

Accept More Element Data Types #22

Open

Firionus added this to the v0.3 milestone Mar 29, 2023

Firionus mentioned this issue Jun 20, 2023

nanmedian might be slower than a naive implementation brenhinkeller/NaNStatistics.jl#41

Closed

Firionus closed this as completed Jun 25, 2023

Firionus mentioned this issue Aug 13, 2023

Update README.md about FastRunningMedian NaN Handling yha/RunningQuantiles.jl#2

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Support for NaN input #21

Support for NaN input #21

wuhuangjianCN commented Mar 22, 2023

Firionus commented Mar 23, 2023 •

edited

Loading

wuhuangjianCN commented Mar 24, 2023 •

edited

Loading

wuhuangjianCN commented Mar 27, 2023

Option 1: Poison Whole Window

Option 2: Use Sorting Convention From Base

Option 3: Use Inverse Sorting Convention From Base

Option 4: Ignore NaN/missing Unless They are the Only Thing in Window

Conclusion

Firionus commented Mar 29, 2023 •

edited

Loading

Firionus commented Mar 29, 2023 •

edited

Loading

wuhuangjianCN commented Mar 30, 2023

berjine commented Jun 19, 2023 •

edited

Loading

Firionus commented Jun 19, 2023

berjine commented Jun 20, 2023 •

edited

Loading

Firionus commented Jun 20, 2023 •

edited

Loading

Firionus commented Jun 20, 2023

berjine commented Jun 20, 2023 •

edited

Loading

Firionus commented Jun 20, 2023

Firionus commented Jun 25, 2023 •

edited

Loading

Support for NaN input #21

Support for NaN input #21

Comments

wuhuangjianCN commented Mar 22, 2023

Firionus commented Mar 23, 2023 • edited Loading

Option 1: Poison Whole Window

Option 2: Use Sorting Convention From Base

Option 3: Use Inverse Sorting Convention From Base

Option 4: Ignore NaN/missing Unless They are the Only Thing in Window

Conclusion

wuhuangjianCN commented Mar 24, 2023 • edited Loading

wuhuangjianCN commented Mar 27, 2023

Option 1: Poison Whole Window

Option 2: Use Sorting Convention From Base

Option 3: Use Inverse Sorting Convention From Base

Option 4: Ignore NaN/missing Unless They are the Only Thing in Window

Conclusion

Firionus commented Mar 29, 2023 • edited Loading

Firionus commented Mar 29, 2023 • edited Loading

wuhuangjianCN commented Mar 30, 2023

berjine commented Jun 19, 2023 • edited Loading

Firionus commented Jun 19, 2023

berjine commented Jun 20, 2023 • edited Loading

Firionus commented Jun 20, 2023 • edited Loading

Firionus commented Jun 20, 2023

berjine commented Jun 20, 2023 • edited Loading

Firionus commented Jun 20, 2023

Firionus commented Jun 25, 2023 • edited Loading

Firionus commented Mar 23, 2023 •

edited

Loading

wuhuangjianCN commented Mar 24, 2023 •

edited

Loading

Firionus commented Mar 29, 2023 •

edited

Loading

Firionus commented Mar 29, 2023 •

edited

Loading

berjine commented Jun 19, 2023 •

edited

Loading

berjine commented Jun 20, 2023 •

edited

Loading

Firionus commented Jun 20, 2023 •

edited

Loading

berjine commented Jun 20, 2023 •

edited

Loading

Firionus commented Jun 25, 2023 •

edited

Loading