mean(img) type instability when img is of eltype Gray{N0f8}

[johnnychen94]

julia> versioninfo()
Julia Version 1.5.0
Commit 96786e22cc (2020-08-01 23:44 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8

julia> using ImageCore, Statistics

julia> x = rand(Gray{N0f8}, 4, 4);

julia> @code_warntype mean(x)

Variables
  #self#::Core.Compiler.Const(Statistics.mean, false)
  A::Array{Gray{Normed{UInt8,8}},2}

Body::Union{}
1 ─     Statistics.:(var"#mean#2")(Statistics.:(:), #self#, A)
└──     Core.Compiler.Const(:(return %1), false)

which is good, but once I load ColorVectorSpace, it isn't anymore:

julia> using ColorVectorSpace

julia> @code_warntype mean(x)
Variables
  #self#::Core.Compiler.Const(Statistics.mean, false)
  A::Array{Gray{Normed{UInt8,8}},2}

Body::Union{Gray{Float32}, Gray{Float64}}
1 ─ %1 = Statistics.:(var"#mean#2")(Statistics.:(:), #self#, A)::Union{Gray{Float32}, Gray{Float64}}
└──      return %1

Both Gray{Float32} and N0f8 are good here.

[kimikage]

Perhaps @timholy already knows about this problem. JuliaMath/FixedPointNumbers.jl#183 is listed in #131 (comment).

Julia v1.5.0 has already been released, so I think it's worth fixing separately.

[timholy]

Thanks for spotting this!

which is good

A return type of Union{} (see Body::Union{}) means it doesn't return:

julia> using ImageCore, Statistics

julia> x = rand(Gray{N0f8}, 4, 4);

julia> mean(x)
ERROR: MethodError: no method matching /(::Gray{Normed{UInt8,8}}, ::Int64)
Closest candidates are:
  /(::Graphics.Vec2, ::Real) at /home/tim/.julia/packages/Graphics/GTY20/src/Graphics.jl:111
  /(::Missing, ::Number) at missing.jl:115
  /(::BigInt, ::Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}) at gmp.jl:540
  ...
Stacktrace:
 [1] _mean(::typeof(identity), ::Array{Gray{Normed{UInt8,8}},2}, ::Colon) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.5/Statistics/src/Statistics.jl:176
 [2] mean(::Array{Gray{Normed{UInt8,8}},2}; dims::Function) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.5/Statistics/src/Statistics.jl:164
 [3] mean(::Array{Gray{Normed{UInt8,8}},2}) at /home/tim/src/julia-1/usr/share/julia/stdlib/v1.5/Statistics/src/Statistics.jl:164
 [4] top-level scope at REPL[3]:1

julia> using ColorVectorSpace

julia> mean(x)
Gray{Float32}(0.5892157f0)

so it's inferrable only as an error.

The possibility of a Gray{Float64} comes from https://github.com/JuliaLang/Statistics.jl/blob/b384104d35ff0e7cf311485607b177223ed72b9a/src/Statistics.jl#L170. Fundamentally, the problem is that we have

julia> x = [0.4N0f8, 0.6N0f8]
2-element Array{N0f8,1} with eltype Normed{UInt8,8}:
 0.4N0f8
 0.6N0f8

julia> sum(x)/2
0.5

julia> sum(x./2)
0.5f0

which seems less than ideal. I suspect the right answer is that we should change the accumulator type to Float32. Thoughts?

[kimikage]

I suspect the right answer is that we should change the accumulator type to Float32.

However, Treduce = Float64 is a 5-year proven specification. (cf. JuliaMath/FixedPointNumbers.jl#31)

Of course, it is possible to define accumulator types for color types separately from FixedPoint, but Gray and Real would be inconsistent.

I'm in favor of rethinking the specifications, but I think it's better to repair the backwards compatibility first. FixedPointNumbers v0.9 will take a bit longer to be released.

[timholy]

However, Treduce = Float64 is a 5-year proven specification. (cf. JuliaMath/FixedPointNumbers.jl#31)

I think your edited comment suggests we pretty much agree on how to proceed, but to be sure there's no misunderstanding let me make sure I've explained my viewpoint. Just because something has been in a code base for a long time does not mean it's right. As we strive for 1.0, we should do what's right. Julia's usage is growing rapidly, and the number of users 10 years from now will dwarf the size of the community now. Since we're not at 1.0, we should aim for as much logical consistency as we can even if we have to make some breaking changes.

[timholy]

But yes, it's fine to make non-breaking changes first, as long as that doesn't start holding back dependent packages in some way. (I can't see how this one would.)

[johnnychen94]

My best guess is that if we consistently use float(T) (or floattype(T)), then this issue just gets fixed automatically.

JuliaMath/FixedPointNumbers.jl#127

[kimikage]

I'm fine with the Float32 accumulator type for 8-/16-bit FixedPoint. However, I don't believe that Float32 is sufficient, as it "only" has 24-bit precision. Also, it is not so desirable for the accumulator type for 64-bit FixedPoint to be BigFloat.

In the first place, I think it is worth doing some specialization in terms of accuracy and speed.
(cf. JuliaMath/FixedPointNumbers.jl#146 (comment)) There is no API in Base to do that safely, though.

Again, this is not only an issue with the nightly build, but also with v1.5.