How far does type instability propagate?

I have a function type_uint(val::Integer) ::Type{<:Unsigned} which returns the smallest UInt type that can still represent val without overflowing. I use this everywhere in my code, but only at top level.
For example, given l\times m and m\times n sparse matrices X and Y, the product will be l\times n, so the indices/positions will be of type tp=type_uint(l) and the values/weights will be of type tw=promote_type(twx,twy). I pass tp(l), tw(0) and the indices and values of X and Y to a function that does all the work:

_compute_product(l::tp, ::tw,
    Xcolptr::Vector{tPx}, Xrowval::Vector{tpx}, Xnzval::Vector{twx}, 
    Ycolptr::Vector{tPy}, Yrowval::Vector{tpy}, Ynzval::Vector{twy}) 
where {tp<:Integer, tw, tPx<:Integer, tpx<:Integer, twx, tPy<:Integer, tpy<:Integer, twy}

Of course, the latter calls a bunch of other functions (which are parametrized with constrained types like above) that run in hot loops, optionally in parallel.

My question: Obviously tp(l)::Any is type unstable. But when I pass it to later functions, its type is a concrete integer. Does type-instability stop at the first function, or does it propagate to other functions? Since later functions do not accept Any arguments, I thought I was safe from losing performance. Is my thinking correct? When I use @code_warntype, I get:

MethodInstance for *(::MatCSC{UInt24, UInt24, Float64}, ::MatCSC{UInt24, UInt24, Float64}, ::Type{mtd3}, ::Type{prl1})
  from *(X::MatCSC{tPx, tpx, twx}, Y::MatCSC{tPy, tpy, twy}, mtd::Type{<:SLATH.Mtd}, prl::Type{<:SLATH.Prl}) where {tPx, tpx, twx, tPy, tpy, twy} @ SLATH ~/venvs_julia/SLATH/src/SLATH.jl:2317
Static Parameters
  tPx = UInt24
  tpx = UInt24
  twx = Float64
  tPy = UInt24
  tpy = UInt24
  twy = Float64
Arguments
  #self#::Core.Const(*)
  X::MatCSC{UInt24, UInt24, Float64}
  Y::MatCSC{UInt24, UInt24, Float64}
  mtd::Core.Const(mtd3)
  prl::Core.Const(prl1)
Locals
  @_6::Int64
  @_7::Int64
  ww::Vector{Float64}
  pp::Vector
  I::Vector{UInt64}
  _0::Float64
  tp::Type
  tP::Type{UInt64}
  n::Int64
  _m::Int64
  m::Int64
  l::Int64
  @_18::MatCSC{UInt64, T, Float64} where T<:Integer
Body::MatCSC{UInt64, T, Float64} where T<:Integer
1 ─        Core.NewvarNode(:(@_6))
β”‚          Core.NewvarNode(:(ww))
β”‚          Core.NewvarNode(:(pp))
β”‚          Core.NewvarNode(:(I))
β”‚          Core.NewvarNode(:(_0))
β”‚          Core.NewvarNode(:(tp))
β”‚          Core.NewvarNode(:(tP))
β”‚   %8   = SLATH.MatCSC::Core.Const(MatCSC)
β”‚   %9   = SLATH.size::Core.Const(size)
β”‚   %10  = (%9)(X)::Tuple{Int64, Int64}
β”‚   %11  = SLATH.size::Core.Const(size)
β”‚   %12  = (%11)(Y)::Tuple{Int64, Int64}
β”‚   %13  = Core._apply_iterate(Base.iterate, Core.tuple, %10, %12)::NTuple{4, Int64}
β”‚   %14  = Base.indexed_iterate(%13, 1)::Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(2)])
β”‚          (l = Core.getfield(%14, 1))
β”‚          (@_7 = Core.getfield(%14, 2))
β”‚   %17  = @_7::Core.Const(2)
β”‚   %18  = Base.indexed_iterate(%13, 2, %17)::Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(3)])
β”‚          (m = Core.getfield(%18, 1))
β”‚          (@_7 = Core.getfield(%18, 2))
β”‚   %21  = @_7::Core.Const(3)
β”‚   %22  = Base.indexed_iterate(%13, 3, %21)::Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(4)])
β”‚          (_m = Core.getfield(%22, 1))
β”‚          (@_7 = Core.getfield(%22, 2))
β”‚   %25  = @_7::Core.Const(4)
β”‚   %26  = Base.indexed_iterate(%13, 4, %25)::Core.PartialStruct(Tuple{Int64, Int64}, Any[Int64, Core.Const(5)])
β”‚          (n = Core.getfield(%26, 1))
β”‚   %28  = SLATH.:(==)::Core.Const(==)
β”‚   %29  = m::Int64
β”‚   %30  = _m::Int64
β”‚   %31  = (%28)(%29, %30)::Bool
└──        goto #3 if not %31
2 ─        goto #4
3 ─ %34  = Base.throw::Core.Const(throw)
β”‚   %35  = Base.AssertionError::Core.Const(AssertionError)
β”‚   %36  = Main.Base::Core.Const(Base)
β”‚   %37  = Base.getproperty(%36, :inferencebarrier)::Any
β”‚   %38  = Main.Base::Core.Const(Base)
β”‚   %39  = Base.getproperty(%38, :string)::Any
β”‚   %40  = (%37)(%39)::Any
β”‚   %41  = m::Int64
β”‚   %42  = _m::Int64
β”‚   %43  = Base.string("X*Y is undefined, since size(X,2) = ", %41, " β‰  ", %42, " = size(Y,1).")::Any
β”‚   %44  = (%40)(%43)::Any
β”‚   %45  = (%35)(%44)::Any
└──        (%34)(%45)
4 β”„        (tP = SLATH.UInt64)
β”‚   %48  = l::Int64
β”‚          (tp = SLATH.type_uint(%48))
β”‚   %50  = SLATH.:*::Core.Const(*)
β”‚   %51  = SLATH.zero::Core.Const(zero)
β”‚   %52  = $(Expr(:static_parameter, 3))::Core.Const(Float64)
β”‚   %53  = (%51)(%52)::Core.Const(0.0)
β”‚   %54  = SLATH.zero::Core.Const(zero)
β”‚   %55  = $(Expr(:static_parameter, 6))::Core.Const(Float64)
β”‚   %56  = (%54)(%55)::Core.Const(0.0)
β”‚          (_0 = (%50)(%53, %56))
β”‚   %58  = SLATH.zeros::Core.Const(zeros)
β”‚   %59  = tP::Core.Const(UInt64)
β”‚   %60  = n::Int64
β”‚   %61  = (%60 + 1)::Int64
β”‚          (I = (%58)(%59, %61))
β”‚   %63  = I::Vector{UInt64}
β”‚          Base.setindex!(%63, 1, 1)
β”‚   %65  = SLATH._mul_CS_CS!::Core.Const(SLATH._mul_CS_CS!)
β”‚   %66  = Base.getproperty(X, :I)::Vector{UInt24}
β”‚   %67  = Base.getproperty(X, :pp)::Vector{UInt24}
β”‚   %68  = Base.getproperty(X, :ww)::Vector{Float64}
β”‚   %69  = Base.getproperty(Y, :I)::Vector{UInt24}
β”‚   %70  = Base.getproperty(Y, :pp)::Vector{UInt24}
β”‚   %71  = Base.getproperty(Y, :ww)::Vector{Float64}
β”‚   %72  = I::Vector{UInt64}
β”‚   %73  = tp::Type
β”‚   %74  = (%73)(0)::Any
β”‚   %75  = _0::Core.Const(0.0)
β”‚   %76  = l::Int64
β”‚   %77  = m::Int64
β”‚   %78  = n::Int64
β”‚   %79  = (%65)(mtd, prl, %66, %67, %68, %69, %70, %71, %72, %74, %75, %76, %77, %78)::Tuple{Vector, Vector{Float64}}
β”‚   %80  = Base.indexed_iterate(%79, 1)::Core.PartialStruct(Tuple{Vector, Int64}, Any[Vector, Core.Const(2)])
β”‚          (pp = Core.getfield(%80, 1))
β”‚          (@_6 = Core.getfield(%80, 2))
β”‚   %83  = @_6::Core.Const(2)
β”‚   %84  = Base.indexed_iterate(%79, 2, %83)::Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(3)])
β”‚          (ww = Core.getfield(%84, 1))
β”‚   %86  = SLATH.MatCSC::Core.Const(MatCSC)
β”‚   %87  = tP::Core.Const(UInt64)
β”‚   %88  = tp::Type
β”‚   %89  = SLATH.typeof::Core.Const(typeof)
β”‚   %90  = _0::Core.Const(0.0)
β”‚   %91  = (%89)(%90)::Core.Const(Float64)
β”‚   %92  = Core.apply_type(%86, %87, %88, %91)::Type{MatCSC{UInt64, T, Float64}} where T
β”‚   %93  = tp::Type
β”‚   %94  = l::Int64
β”‚   %95  = Base.getindex(%93, %94)::Vector
β”‚   %96  = I::Vector{UInt64}
β”‚   %97  = pp::Vector
β”‚   %98  = ww::Vector{Float64}
β”‚   %99  = (%92)(%95, %96, %97, %98)::MatCSC{UInt64, T, Float64} where T<:Integer
β”‚          (@_18 = %99)
β”‚   %101 = @_18::MatCSC{UInt64, T, Float64} where T<:Integer
β”‚   %102 = (%101 isa %8)::Core.Const(true)
└──        goto #6 if not %102
5 ─        goto #7
6 ─        Core.Const(:(@_18))
β”‚          Core.Const(:(Base.convert(%8, %105)))
└──        Core.Const(:(@_18 = Core.typeassert(%106, %8)))
7 β”„ %108 = @_18::MatCSC{UInt64, T, Float64} where T<:Integer
└──        return %108

Is this healthy?

Off-topic, a second question that I have, how can I force julia to specialize in the following case:

function Base.hcat(XX::MatCSC...)  @assert length(XX)>=3; ...general implementation ... end
function Base.hcat(X::MatCSC{tP,tp,tw}) where {...}  return X end
function Base.hcat(X::MatCSC{tPx,tpx,twx}, Y::MatCSC{tPy,tpy,twy}) where {...}  ...more efficient implementation ... end;

When calling on 1 or 2 arguments, I get the assertion error, even though the 1- and 2- case is defined later and is more specific (parametrized by types).

Yes, this is the notion of function barriers.

Are you sure the type parameters are correct in the 1 or 2 argument call?

3 Likes

Ok, good, my later functions are safe from type instability. What about my stack trace of @code_warntype?

Why would type parameters in the 1 or 2 argument call not be correct? I’m calling the hcat method with 1, 2, or multi arguments, all of which are MatCSC.

It’s a bit hard to say without the source code, but there do seem to be a lot of Anys. Most of them look related to some error message, so that’s probably fine? Now,

β”‚   %73  = tp::Type
β”‚   %74  = (%73)(0)::Any

seems a bit weird, though. In any case, you could also just see whether the performance is worse and/or you get more allocations than you would expect.

Because you shouldn’t have to force Julia to specialise in this case, it does it automatically:

julia> f(::Vector...) = println("General")
f (generic function with 1 method)

julia> f(::Vector{T}) where T = println("Specialised 1")
f (generic function with 2 methods)

julia> f(::Vector{S}, ::Vector{T}) where {S, T} = println("Specialised 2")
f (generic function with 3 methods)

julia> f(rand(1))
Specialised 1

julia> f(rand(1), rand(2))
Specialised 2

julia> f(rand(Int32, 1), rand(2))
Specialised 2

julia> f(rand(1), rand(2), rand(3))
General

So my best guess is that e.g. in your single argument assertion error case, you’re simply not supplying a MatCSC{tP,tp,tw}.

Your example with f was helpful. I figured out that by changing
function Base.hcat(XX::MatCSC...) ::MatCSC to
function Base.hcat(XX::MatCSC{<:Integer,<:Integer,Any}...) ::MatCSC, so that all 3 methods are parametrized, causes the right specialization. But why does this happen in my hcat and not in your f?

Even more strangely, using <:Any instead of Any raises the assertion error again.

Could you provide a MWE?

Those functions should not be more specialized, MatCSC and MatCSC{...} where {...} should actually be equivalent as stated in the docs More about types Β· The Julia Language

For example

julia> Array === (Array{T} where T)
true

I can’t see why can’t you use a single implementation with varaargs and then call a secondary function

function _foo1(x) ... end
function _foo2(x1, x2) ... end
function foo(xx...)
    length(xx) == 1 && _foo1(xx[1])
    length(xx) == 2 && _foo2(xx[1], xx[2])
    ... your implementation...
end

This should not create performance problems. If somehow you still want to specilize the functions, a more β€œhacky” way could be to use NTuples

julia> foo(x::NTuple{1, Float64}) = println("1 $x")
foo (generic function with 1 method)

julia> foo(x::NTuple{2, Float64}) = println("2 $x")
foo (generic function with 2 methods)

julia> foo(x::NTuple{N, Float64}) where N= println("generic $x")
foo (generic function with 3 methods)

julia> foo(x::Float64...) = foo(NTuple(x))
foo (generic function with 4 methods)

julia> foo(1.)
1 (1.0,)

julia> foo(1., 2.)
2 (1.0, 2.0)

julia> foo(1., 2., 3.)
generic (1.0, 2.0, 3.0)
1 Like

Could you provide a MWE?

Those functions should not be more specialized, MatCSC and MatCSC{…} where {…} should actually be equivalent as stated in the docs

Here is the MWE:

using SparseArrays
struct MatCSC{tP<:Integer, tp<:Integer, tw}
   d ::Vector{tp}
   I ::Vector{tP}
   pp::Vector{tp}
   ww::Vector{tw} 
end
function matCSC(X::SparseMatrixCSC{tw,tp}) ::MatCSC  where {tp<:Integer, tw}
   return MatCSC([size(X,1)], X.colptr, X.rowval, X.nzval) 
end
function f(X::MatCSC{tP,tp,tw})  where {tP,tp,tw}  return X end
function g(X::MatCSC{tP,tp,tw})  where {tP,tp,tw}  return X end
function h(X::MatCSC{tP,tp,tw})  where {tP,tp,tw}  return X end
function f(XX::MatCSC...)                             @assert length(XX) β‰₯ 2;  return XX[1] end
function g(XX::MatCSC{<:Integer,<:Integer,Any}...)    @assert length(XX) β‰₯ 2;  return XX[1] end
function h(XX::MatCSC{<:Integer,<:Integer,<:Any}...)  @assert length(XX) β‰₯ 2;  return XX[1] end
x=sprand(5,5,0.4); X=matCSC(x); 

julia> f(X)
ERROR: AssertionError: length(XX) β‰₯ 2
Stacktrace:
 [1] f(XX::MatCSC{Int64, Int64, Float64})
   @ Main ./REPL[7]:1
 [2] top-level scope
   @ REPL[11]:1

julia> g(X)
MatCSC{Int64, Int64, Float64}([5], [1, 4, 7, 9, 10, 12], [1, 2, 5, 2, 3, 4, 1, 3, 2, 2, 5], [0.21296663444239972, 0.42541787693734945, 0.21276750087826435, 0.3427309883352708, 0.1808337277144113, 0.6450865248125436, 0.9403180620312582, 0.30062699286008066, 0.19106133347617826, 0.7349450918800273, 0.7130889414916233])

julia> h(X)
ERROR: AssertionError: length(XX) β‰₯ 2
Stacktrace:
 [1] h(XX::MatCSC{Int64, Int64, Float64})
   @ Main ./REPL[9]:1
 [2] top-level scope
   @ REPL[13]:1

It seems MatCSC and MatCSC{tP,tp,tw} are not equivalent as arguments to functions. The above example is weird.

One is a strict subtype of the other.

julia> (MatCSC{a,b,c} where {a,b,c}) == MatCSC
false

julia> (MatCSC{a,b,c} where {a,b,c}) <: MatCSC
false

julia> MatCSC <: (MatCSC{a,b,c} where {a,b,c})
true

That’s because your definition itself had type bounds, which just don’t print out fully (anymore). It’s still part of the type, which we can print out with some internals trickery to identify a fully written, equivalent type:

julia> MatCSC.body.body.body
MatCSC{tP<:Integer, tp<:Integer, tw}

julia> MatCSC{<:Integer, <:Integer, <:Any} == MatCSC
true

Now that we’re reminded of what MatCSC really means, the dispatches can be explained:

  1. f(X) dispatches to the now plainly more specific Vararg method, which does accept 1 argument.
  2. g(X) dispatches to the parametric method simply because X is NOT an instance of MatCSC{<:Integer,<:Integer,Any}. The omitted <: turning a type bound into a parameter value is a huge difference, the same way Vector{Float64} is a subtype of Vector{<:Any}, not Vector{Any} (type parameters are invariant).
  1. h(X) dispatches to the Vararg method for the same reasons that f(X) did; the only difference is you wrote out the definition’s type bounds. If you check the methods, you will find that those type bounds aren’t printed out either:
julia> methods(h)
# 2 methods for generic function "h" from Main:
 [1] h(XX::MatCSC...)
     @ REPL[6]:1
 [2] h(X::MatCSC{tP, tp, tw}) where {tP, tp, tw}
     @ REPL[9]:1

It’s good practice to explicitly repeat definitions’ type bounds in subtypes and annotations when you need to specify parameters but don’t have more specific bounds. You might be able to get away with defaulting to <:Any sometimes, but that evidently has its limits.

3 Likes

That MWE shows a behavior but does not really explain why you want/need to do that. If it is a matter of algorithms the best way would be to create a custom type for each algorithm supported and specialize functions based on the algorithm.

The difference in behavior has been well explained by @Benny, and it is given by the type bounds, It can also be seen by calling dump on the type, which explicitly report the bounds

julia> dump(MatCSC{tP,tp,tw} where {tP,tp,tw})
UnionAll
  var: TypeVar
    name: Symbol tP
    lb: Union{}
    ub: Any
  body: UnionAll
    var: TypeVar
      name: Symbol tp
      lb: Union{}
      ub: Any
    body: UnionAll
      var: TypeVar
        name: Symbol tw
        lb: Union{}
        ub: Any
      body: MatCSC{tP, tp, tw} <: Any
        d::Vector{tp}
        I::Vector{tP}
        pp::Vector{tp}
        ww::Vector{tw}

julia> dump(MatCSC{<:Integer,<:Integer})
UnionAll
  var: TypeVar
    name: Symbol #s15
    lb: Union{}
    ub: Integer <: Real
  body: UnionAll
    var: TypeVar
      name: Symbol #s16
      lb: Union{}
      ub: Integer <: Real
    body: UnionAll
      var: TypeVar
        name: Symbol tw
        lb: Union{}
        ub: Any
      body: MatCSC{var"#s15"<:Integer, var"#s16"<:Integer, tw} <: Any
        d::Vector{var"#s16"}
        I::Vector{var"#s15"}
        pp::Vector{var"#s16"}
        ww::Vector{tw}
1 Like

Relevant discussion on Github:

1 Like