Tuple indexing taking time?

A fairly general question about a function that I’m working on before I try to summarize what it’s doing. The profile view is below. There are three main red blocks that take up most of the time, each for the same line of code, which makes sense. But, I don’t understand why there is so much tuple getindex time above these blocks or why there is also empty time above them. Is this telling me that I have a dynamic portion that needs to be fixed? The rest of the operations aren’t trivial, so I don’t have any reason to expect tuple indexing to be a slow part.

Zooming in, the little spike above the getindex block is PyCall.jl, pydecref, which I don’t know how to interpret.

Are you using indexing with multiple values? E.g., t[2:4]? It’d be helpful if you could show a bit of your code.

1 Like

The main red blocks point to the inside of the for loop in this function, where σ += gas[i](T, P, idx):

function _schwarzschild(P::Real, I::Real, param)::Float64
    #unpack
    gas, idx, ν, g, μ, Tfun = param
    #unravel the log pressure coordinate
    P = exp(-P)
    #compute temperature from given function
    T = Tfun(P)
    #get concentration-scaled cross-section from interpolators
    σ = 0.0
    for i = 1:length(gas)
        σ += gas[i](T, P, idx)
    end
    #compute dI/dlogP
    P*schwarzschild(I, ν, σ, g, μ, T)
end

the gas variable is a Tuple but I don’t know why indexing it would take roughly as much time as the operation that each element performs, which is a 2D interpolation. The gas tuple does have elements with different types, though.

Reals are abstract. Not good. Better:

function _schwarzschild(P, I, param)::Float64
1 Like

Does that matter though? My “P::Real, I::Real” in the definition is just supposed to constrain input types, even though P and T are always Float64 in this case.

Also: unpacking param may result in type instability, I think. Function barrier?

It’s argument annotation, so it is not very important (just some modest constraints on argument types).

I guess this tuple is a tuple of functions, and probably program take some run time to figure out, which function to run.

@code_warntype can help here.

Ok I’ll learn how to use @code_warntype and report back.

No, the use of ::Real in a function argument type has no performance penalty whatsoever. Writing function _schwarzschild(P::Real, I::Real, ...). is perfectly fine.

4 Likes

Right you are.

I think what tripped me up was this observation of mine from about a week ago:

julia>  f(2.0)
-2.18504e+00

julia>  f(2)
-2.18504e+00

julia> methods(f)
# 1 method for generic function "f":
[1] f(x::Real) in Main at REPL[5]:1

I noticed that specialized methods had not been derived for f. How do we explain that?

That’s not what methods shows anyway–it just tells you that there is a single method for f, and that method accepts any subtype of ::Real. The actual specialization for a particular type happens when you call f, and you can see that in the @code_warntype f(2)

2 Likes

Didn’t it used to be that I could see those methods when the function was called for different arguments? I think I put an example like this into my tutorial a while back.
Edit: Now I remember, I have to type the argument explicitly.

Not as far as I know (and I just checked that it doesn’t happen in versions as old as Julia 0.5). Maybe you actually defined foo(x::Real) and foo(x::Int) separately in your tutorial?

3 Likes

So, one solution seems to be replacing that gas tuple with a new type. The tuple has different types of objects in it. If I replace the tuple with a single type that explicitly organizes the different gasses, then use a new function on it instead of the for loop, the getindex blocks go away.

I called the new type GasPack. There is still a little bit of empty space above that block though :woman_shrugging:.

1 Like

So would it be correct to say that f(x::Real) = ... was equivalent to f(x::T) where {T<:Real} = ...?

1 Like

I don’t really understand how to decipher the @code_warntype results though. This is the first bit

julia> @code_warntype _schwarzschild(1e5, 1, ((CO2,H2O), 10, 667.0, 9.8, 0.03, M))
Variables
  #self#::Core.Compiler.Const(_schwarzschild, false)
  P@_2::Float64
  I::Int64
  param::Tuple{Tuple{MinorGas,VariableGas{var"#1#2"}},Int64,Float64,Float64,Float64,MoistAdiabat}
  @_5::Int64
  gas::Tuple{MinorGas,VariableGas{var"#1#2"}}
  idx::Int64
  ν::Float64
  g::Float64
  μ::Float64
  Tfun::MoistAdiabat
  T::Float64
  σ::Float64
  @_14::Union{Nothing, Tuple{Int64,Int64}}
  i::Int64
  P@_16::Float64

what is @_14?

I think that’s a temp variable found somewhere in the body below. Could it be that MinorGas is not a completely concrete type? A tuple should work as well as a struct with two entries in terms of compilation and type inference, but if your gas types are abstract or the function you run when calling them is not completely specified by the type, then that’s a type instability

I think this is the basic problem. Indexing into a tuple with different types is generally not type-stable because the type of the result depends on the value of i.

I wonder if GitHub - cstjean/Unrolled.jl: Unrolling loops at compile-time would help here?

5 Likes

That seems right.

It also speeds up if, instead of a tuple, the gas argument is just a Vector{AbstractGas}, where I’ve defined the function for AbstractGas that is called inside the for loop of _schwarzschild.

I think so. They only stop being equivalent when you have them inside type parameters (which are invariant):

julia> struct Wrapper{T}; a :: T; end

julia> f(x :: Wrapper{Real}) = x.a
f (generic function with 1 method)

julia> g(x :: Wrapper{T}) where {T<:Real} = x.a
g (generic function with 1 method)

julia> x1 = Wrapper{Real}(1.0)
Wrapper{Real}(1.0)

julia> x2 = Wrapper{Float64}(1.0)
Wrapper{Float64}(1.0)

julia> f(x1)
1.0

julia> f(x2)
ERROR: MethodError: no method matching f(::Wrapper{Float64})
Closest candidates are:
  f(::Wrapper{Real}) at REPL[2]:1
Stacktrace:
 [1] top-level scope at REPL[7]:1

julia> g(x1)
1.0

julia> g(x2)
1.0

In this case the short way to write g would be:

julia> h(x :: Wrapper{<:Real}) = x.a
h (generic function with 1 method)

julia> h(x1)
1.0

julia> h(x2)
1.0