What crazy instability is this?

joa-quim · May 28, 2020, 2:13pm

The MWE

function img16to8(mat::Array{UInt16}; kw...)
	d = Dict{Symbol,Any}(kw)
    img = Array{UInt8}(undef,size(mat));
    if (haskey(d, :histo_bounds))
        val_ = d[:histo_bounds]
        val = UInt16.(val_)		# This also tests for good bounding values
        val = Int(val)          # Make it same type as when setting it manually below
		len = length(val_)
		if (len == 1)
			#val = 0
			sc = 255 / (65535 - val)
			#@show(typeof(val), typeof(sc))
			@inbounds for k = 1:length(img)
				img[k] = (mat[k] < val) ? 0 : round(UInt8, (mat[k] - val) * sc)
			end
		end
	end
	return img
end

Very slow

img16 = rand(UInt16, 7000, 7000);

@btime img16to8(img16, histo_bounds=0);
  7.348 s (243087474 allocations: 3.67 GiB)

Now, if I remove the comment in #val = 0 in the function, which will do exactly the same operation, and reload it

julia> @btime img16to8(img16, histo_bounds=0);
  82.211 ms (11 allocations: 46.73 MiB)

???

There is clearly a type instability or global variable somewhere but can’t find it. Uncommenting the @show line shows that val & sc have the same type in both tests.

stevengj · May 28, 2020, 2:25pm

In general, you should avoid changing the type of a variable within a function. Sometimes the compiler can optimize it away, but other times not, which is presumably what you are seeing here. Also note that this code will fail if val_ is an array, because Int(...) doesn’t work for arrays.

Maybe you want something like:

vals = d[:histo_bounds]
if length(vals) == 1
    val = Int(first(vals))
    ...
end

joa-quim · May 28, 2020, 2:33pm

Thanks. The val = Int(val) is in this example code only to show that the types of val and sc are the same whether I use the fast or the slow version. It’s not in the original code. But it doesn’t make a difference either.
The val = UInt16.(val_) was an attempt to make the operations in

img[k] = (mat[k] < val) ? 0 : round(UInt8, (mat[k] - val) * sc)

of the same type as closest. I have tried dozens of variations but the huge difference in run times remains.

mbauman · May 28, 2020, 2:37pm

Here’s another: this is sometimes an Int and sometimes a UInt8. Use:

(mat[k] < val) ? UInt8(0) : round(UInt8, (mat[k] - val) * sc)

joa-quim · May 28, 2020, 2:42pm

Had tried that too. But nope, same result.

benninkrs · May 28, 2020, 2:50pm

Confirming that the types are the same at run time is not sufficient. What matters is whether the compiler can determine the types at compile time. Have you tried using @code_warntype ?

mbauman · May 28, 2020, 2:51pm

Your root issue is almost certainly stemming from the kwargs and the Any dict you put them in. Julia won’t specialize on either — and thus it’ll always be ::Any unless you change around what val is. Resetting val=0 will “recover” that instability because from there on out Julia will know that it’s an Int. Doing a strange song and dance hopping between UInt16 and Int may be complicated enough that Julia doesn’t figure it out.

joa-quim · May 28, 2020, 3:05pm

I had actually tried around that line (the dict with Any is something I need to use in the real case) I tried things like

val_ = d[:histo_bounds]
val = parse(Float64, @sprintf("%d", val_))

with mixed results. I mean sometimes it seemed to work but later not.
At this time is works but seems a bit ugly.

mbauman · May 28, 2020, 3:10pm

This is the same idea but is much much simpler. Julia won’t know what vals is, but it’ll know that val must be an Int.

joa-quim · May 28, 2020, 3:21pm

… but it doesn’t improve. Still very slow.

pixel27 · May 28, 2020, 7:04pm

Ultimately it appears to come down to the Dict{Symbol, Any}. That appears to be throwing the compiler for a loop. I’ve reduced your MWE to this (it has less distractions):

function img16to8(mat::Array{UInt16}; kw...)
	d   = Dict{Symbol, Any}(kw)
    img = zeros(UInt8, size(mat))
	val = UInt16(d[:histo_bounds])
	sc  = 255 / (65535 - val)
	@inbounds for k = 1:length(img)
		if mat[k] >= val
			tmp = (mat[k] - val) * sc
			img[k] = round(UInt8, tmp)
		end
	end
	return img
end

img16to8(rand(UInt16, 7000, 7000); histo_bounds=0)

Running this with julia --track-allocation=all I get the results:

        - function img16to8(mat::Array{UInt16}; kw...)
  4333364 	d   = Dict{Symbol, Any}(kw)
 49000144   img = zeros(UInt8, size(mat))
        0 	val = UInt16(d[:histo_bounds])
       32 	sc  = 255 / (65535 - val)
        0 	@inbounds for k = 1:length(img)
        0 		if mat[k] >= val
1555774112 			tmp = (mat[k] - val) * sc
        0 			img[k] = round(UInt8, tmp)
        - 		end
        - 	end
        0 	return img
        - end
        - 
        - img16to8(rand(UInt16, 7000, 7000); histo_bounds=0)
        -

The allocations at line 8 (tmp = (mat[k] - val) * sc) is what is causing your slow down. Julia doesn’t appear to be doing a “normal” calculation, but treating val as an any with all the overhead that entails. If you change the Dict to Dict{Symbol, Int)(kw) the compiler generates code that doesn’t create the allocations.

I would have thought the forcing val to be a UInt16 would allow the compiler to generate efficient code but no such luck.

pixel27 · May 28, 2020, 7:17pm

If I go to this, which still exhibits the issue, but seems to show a better result with @code_warntype:

function img16to8(mat::Array{UInt16}; histo_bounds = 0)
	d   = Dict{Symbol, Any}(:histo_bounds => histo_bounds)
    img = zeros(UInt8, size(mat))
	val = convert(UInt16, d[:histo_bounds])
	sc  = 255 / (65535 - val)
	@inbounds for k = 1:length(mat)
		if mat[k] >= val
			tmp = (mat[k] - val) * sc
			img[k] = round(UInt8, tmp)
		end
	end
	return img
end

@code_warntype img16to8(rand(UInt16, 7000, 7000); histo_bounds=0)

The output doesn’t appear to indicate any type instability:

  #unused#::Core.Compiler.Const(var"#img16to8##kw"(), false)
  @_2::NamedTuple{(:histo_bounds,),Tuple{Int64}}
  @_3::Core.Compiler.Const(img16to8, false)
  mat::Array{UInt16,2}
  histo_bounds::Int64
  @_6::Int64

Body::Array{UInt8,2}
1 ─ %1  = Base.haskey(@_2, :histo_bounds)::Core.Compiler.Const(true, false)
│         %1
│         (@_6 = Base.getindex(@_2, :histo_bounds))
└──       goto #3
2 ─       Core.Compiler.Const(:(@_6 = 0), false)
3 ┄       (histo_bounds = @_6)
│   %7  = (:histo_bounds,)::Core.Compiler.Const((:histo_bounds,), false)
│   %8  = Core.apply_type(Core.NamedTuple, %7)::Core.Compiler.Const(NamedTuple{(:histo_bounds,),T} where T<:Tuple, false)
│   %9  = Base.structdiff(@_2, %8)::Core.Compiler.Const(NamedTuple(), false)
│   %10 = Base.pairs(%9)::Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false)
│   %11 = Base.isempty(%10)::Core.Compiler.Const(true, false)
│         %11
└──       goto #5
4 ─       Core.Compiler.Const(:(Base.kwerr(@_2, @_3, mat)), false)
5 ┄ %15 = Main.:(var"#img16to8#11")(histo_bounds, @_3, mat)::Array{UInt8,2}
└──       return %15

joa-quim · May 28, 2020, 7:30pm

Thanks, I’m using a Any because this chunk came from a larger code where I need it to be a Any. But for this articular case I can live well setting it to Int … and at same time solving the issue.

tim.holy · May 28, 2020, 11:04pm

You’re not quite forcing it:

julia> methods(convert, (Type{UInt16}, Any))
# 7 methods for generic function "convert":
[1] convert(::Type{T}, x::T) where T<:Number in Base at number.jl:6
[2] convert(::Type{T}, x::Number) where T<:Number in Base at number.jl:7
[3] convert(::Type{T}, x::Ptr) where T<:Integer in Base at pointer.jl:23
[4] convert(::Type{T}, x::Base.TwicePrecision) where T<:Number in Base at twiceprecision.jl:250
[5] convert(::Type{T}, x::AbstractChar) where T<:Number in Base at char.jl:180
[6] convert(::Type{T}, index::CartesianIndex{1}) where T<:Number in Base.IteratorsMD at multidimensional.jl:135
[7] convert(::Type{T}, x::T) where T in Base at essentials.jl:171

One of these 7 infers with a return value of Any, so the compiler has to conservatively choose Any. You can write val = convert(UInt16, d[:histo_bounds])::UInt16 to declare that the conversion to UInt16 will actually succeed.

convert is just an ordinary method, there is no language-level guarantee that it returns an object of the type you’ve requested.

Palli · May 28, 2020, 11:22pm

Maybe you got your answer (I didn’t read rest of the thread).

FYI, to debug such (I’m trying for the first time):

julia> @time using Traceur
[ Info: Precompiling Traceur [37b6cedf-1f77-55f8-9503-c64b63398394]
  2.970901 seconds (1.24 M allocations: 64.574 MiB)

julia> @trace img16to8(img16, histo_bounds=0);
┌ Warning: x is assigned as Union{Nothing, Tuple{Symbol,Int64}}
└ @ iterators.jl:234

 [seemingly hung here, then kept on going]

┌ Warning: val is assigned as Union{Nothing, Dict{Symbol,Any}}
└ @ dict.jl:388
┌ Warning: x is assigned as Union{Nothing, Tuple{Symbol,Int64}}
└ @ iterators.jl:234

[something missing before "is"...?]

┌ Warning:  is assigned as Union{Nothing, Tuple{Pair{Symbol,Int64},Int64}}
└ @ dict.jl:102
┌ Warning:  is assigned as Union{Nothing, Tuple{Pair{Symbol,Int64},Int64}}
└ @ dict.jl:103
┌ Warning: getindex returns Any
└ @ array.jl:799
┌ Warning: getindex returns Any
└ @ dict.jl:465

[still waiting here]

joa-quim · May 29, 2020, 1:16am

Thanks, I’ve used traceur sometimes but not easy to interpret all its output.

I have the issue basically solved, but the real case still shows the instability time to time. Then after some restarts it goes away.

pixel27 · May 29, 2020, 2:22am

From the manual: Conversion and Promotion · The Julia Language

The convert function generally takes two arguments: the first is a type object and the second is a value to convert to that type. The returned value is the value converted to an instance of given type.

And then:

Conversion isn’t always possible, in which case a MethodError is thrown indicating that convert doesn’t know how to perform the requested conversion:

So I just sort of thought convert would convert to the object type I told it to and throw an exception if it could not. I’m not quite sure how to interpret what happens based on what you are telling me…

tim.holy · May 29, 2020, 6:40am

The manual and docstring describes the interface it should have, but remember most Julia functions are extensible. For example, I can do this:

julia> struct EmptyType end

julia> Base.UInt16(x::EmptyType) = "hello"

julia> Base.convert(::Type{UInt16}, x::EmptyType) = 1.0

julia> x = EmptyType()
EmptyType()

julia> UInt16(x)
"hello"

julia> convert(UInt16, x)
1.0

That’s terrible design, but nothing at the language level prevents it. (A frequently-requested feature sometimes called “interfaces” or “protocols” would make it possible to ban such methods.)

Let me also correct something I implied that’s a bit wrong. I implied that UInt16(x) called convert(UInt16, x), but actually it now the other way around (it used to be the opposite, and I’m slow to change…). So the UInt16(val::Any) calls one of these:

julia> methods(UInt16, (Any,))
# 13 methods for type constructor:
[1] UInt16(x::Union{Bool, Int32, Int64, UInt32, UInt64, UInt8, Int128, Int16, Int8, UInt128, UInt16}) in Core at boot.jl:711
[2] UInt16(x::Float32) in Base at float.jl:685
[3] UInt16(x::Float64) in Base at float.jl:685
[4] (::Type{T})(x::Float16) where T<:Integer in Base at float.jl:71
[5] (::Type{T})(z::Complex) where T<:Real in Base at complex.jl:37
[6] (::Type{T})(x::Rational) where T<:Integer in Base at rational.jl:109
[7] (::Type{T})(x::BigInt) where T<:Union{UInt128, UInt16, UInt32, UInt64, UInt8} in Base.GMP at gmp.jl:347
[8] (::Type{T})(x::BigFloat) where T<:Integer in Base.MPFR at mpfr.jl:324
[9] (::Type{T})(x::T) where T<:Number in Core at boot.jl:716
[10] (::Type{T})(x::Base.TwicePrecision) where T<:Number in Base at twiceprecision.jl:243
[11] (::Type{T})(x::AbstractChar) where T<:Union{AbstractChar, Number} in Base at char.jl:50
[12] (::Type{T})(x::Enum{T2}) where {T<:Integer, T2<:Integer} in Base.Enums at Enums.jl:19
[13] (dt::Type{var"#s827"} where var"#s827"<:Integer)(ip::Sockets.IPAddr) in Sockets at /home/tim/src/julia-master/usr/share/julia/stdlib/v1.6/Sockets/src/IPAddr.jl:11

First, you’ll note there is no specific method for x::Any. Moreover, in inference Any means “I don’t know” and not “it’s of type Any” since no actual object has type Any. Consequently inference has to deal with the possibility that any of a whole list of methods might be called. For reasons of performance, Julia does not run inference on a method/type combination until necessary, so Julia doesn’t know in advance that all 13 of those methods either return a UInt16 or throw an error. If the number of methods is 4 or fewer, Julia will run inference on all of them and might discover that they all return an UInt16, but beyond 4 methods it instead sets the return type as Any. The tradeoffs are discussed in new call site inference algorithm · Issue #34742 · JuliaLang/julia · GitHub.

lwabeke · May 29, 2020, 9:31am

@pixel27 version on my machine:

julia> @time img16to8(rand(UInt16, 7000, 7000); histo_bounds=0);
  7.978230 seconds (242.70 M allocations: 3.753 GiB, 6.99% gc time)

You can fix this by forcing the type of val on line 4:

julia> function img16to8_quick(mat::Array{UInt16}; kw...)
               d   = Dict{Symbol, Any}(kw)
               img = zeros(UInt8, size(mat))
               val :: UInt16 = UInt16(d[:histo_bounds])
               sc  = 255 / (65535 - val)
               @inbounds for k = 1:length(img)
                       if mat[k] >= val
                               tmp = (mat[k] - val) * sc
                               img[k] = round(UInt8, tmp)
                       end
               end
               return img
       end
img16to8_quick (generic function with 1 method)

julia> @time img16to8_quick(rand(UInt16, 7000, 7000); histo_bounds=0);
  0.212195 seconds (36.16 k allocations: 142.041 MiB, 9.58% gc time)

julia> @time img16to8_quick(rand(UInt16, 7000, 7000); histo_bounds=0);
  0.204825 seconds (20 allocations: 140.191 MiB)

You can even remove the call to UInt16(), it seems to give the same result:

val :: UInt16 = (d[:histo_bounds])

joa-quim · May 29, 2020, 2:31pm

Thanks to all that contributed to enlightening the tricky situations on converting types when a Any is involved.