fix implicit conversion of range to array

[MilesCranmer]

std is inadvertently rerouting std(::AbstractRange) through var(::AbstractArray) rather than var(::AbstractRange), since the latter does not support keywords, and the bug JuliaLang/julia#9498 exists (patched in JuliaLang/julia#60499, which led me to spot this downstream effect)

For example:

julia> @b 100_000_000 std(1:_)
91.392 ms

julia> @b 100_000_000 var(1:_)
5.610 ns

This PR fixes this so that std(::AbstractRange) now uses var(::AbstractRange)

[codecov-commenter]

Codecov Report

❌ Patch coverage is 83.33333% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 96.43%. Comparing base (22dee82) to head (2556a94).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Statistics.jl	83.33%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #197      +/-   ##
==========================================
- Coverage   96.66%   96.43%   -0.24%     
==========================================
  Files           2        2              
  Lines         450      449       -1     
==========================================
- Hits          435      433       -2     
- Misses         15       16       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[MilesCranmer]

@nalimilan would you mind reviewing this?

[nalimilan]

Good catch!

Are you sure existing tests cover all cases in the new code? I see one line which isn't covered in Codecov and I wonder whether there might be corner cases.

[MilesCranmer]

Semantically, we shouldn't have to call abs2 for a method which doesn't take complex numbers (since it's a range).

I just looked this up and actually AbstractRange doesn't constrain to reals:

julia> LinRange(0.0im, 1.0im + 1.0, 10)
10-element LinRange{ComplexF64, Int64}:
 0.0+0.0im, 0.111111+0.111111im, 0.222222+0.222222im, …, 0.888889+0.888889im, 1.0+1.0im

julia> typeof(LinRange(0.0im, 1.0im + 1.0, 10)) <: AbstractRange
true

I'll switch back to abs2

[MilesCranmer]

Let me know if you'd prefer this:

function _var(A::AbstractArray, corrected::Bool, mean, dims::Dims) where Dims
    if mean === nothing
        mean = dims === (:) ? Statistics.mean(A) : Statistics.mean(A, dims=dims)
    end
    v = if dims === (:)
        varm(A, mean; corrected=corrected)
    else
        varm(A, mean; corrected=corrected, dims=dims)
    end
    return dims === (:) ? real(v) : v
end

or keep it split as-is

[nalimilan]

What I would do is this:

function _var(A::AbstractArray, corrected::Bool, mean, dims::Colon)
    if mean === nothing
        mean = Statistics.mean(A)
    end
    return varm(A, mean; corrected=corrected)
end

function _var(A::AbstractArray, corrected::Bool, mean, dims::Dims) where Dims
    if mean === nothing
        mean = Statistics.mean(A, dims=dims)
    end
    return real(varm(A, mean; corrected=corrected, dims=dims))
end

Even if that's a bit verbose, the advantage is that it's clear that it calls varm, so no need to check whether the methods give the same results (which is always hard in corner cases).

[nalimilan]

Why change this? The current method has the advantage that we always return a value of the type of vv, and it's simpler.

[MilesCranmer]

It doesn't change the behavior; just avoids redundant calculations

[nalimilan]

I mean, why not keep return typeof(vv)(NaN)? We don't care about redundant calculations in that case. Also, I suspect the compiler is able to eliminate any of these.

[MilesCranmer]

vv was moved to underneath this branch now to avoid the redundant calculation, since all we want here is the return type. I can revert it if you want though

[nalimilan]

Yeah we don't care about performance in this corner case, and the compiler will probably be able to eliminate all of these.

[nalimilan]

Mmm but maybe we don't even need these functions, thanks to the already existing _var(A::AbstractArray, corrected::Bool, mean, dims) and _var(A::AbstractArray, corrected::Bool, mean, ::Colon) which do exactly that?

[MilesCranmer]

Yeah we've gone full circle; this is the current code, untouched by this PR

function _var(A::AbstractArray, corrected::Bool, mean, dims)
  if mean === nothing
      mean = Statistics.mean(A, dims=dims)
  end
  return varm(A, mean; corrected=corrected, dims=dims)
end

function _var(A::AbstractArray, corrected::Bool, mean, ::Colon)
  if mean === nothing
      mean = Statistics.mean(A)
  end
  return real(varm(A, mean; corrected=corrected))
end

Adding ::Dims to the first one would be equivalent specialization

[nalimilan]

Sorry at first I thought the methods you proposed used ::AbstractRange, not ::AbstractArray. Then can't we just keep them as they are, and just drop var(v::AbstractRange)?

[MilesCranmer]

I am not sure I understand

[nalimilan]

I mean just remove var(v::AbstractRange) and rely on the existence of the _var(A::AbstractArray, ... methods so that var on a range ends up calling varm.

[MilesCranmer]

Maybe this refactoring could be done in the future? I only wanted to fix the one issue I happened upon, I hadn't planned to spend this much time on it.

[nalimilan]

But that's not a refactoring, it's just removing a function. Or do you mean that doesn't work? What you've done on var(v::AbstractRange) seems actually more work than just removing it.

[MilesCranmer]

Well if we fully remove the var(v::AbstractRange) branch then we lose the fast formula abs2(step(v)) * (l+1)*l/12 so I'm not quite sure I follow

[nalimilan]

But we can still use if m == nothing (either with an if or with a separate method m::Nothing). Really, varm is just var(... mean=...) designed before keyword arguments were fast, the two are exactly the same.

[MilesCranmer]

Fixed!

[MilesCranmer]

Implemented

[MilesCranmer]

Ok thanks updated

[nalimilan]

Thanks! So now we get more features almost without more code!