Avoiding allocations for function arguments with non-concrete types

[Nicholaswogan]

Hi, I am beginning to learn how to build my own types, and pass them to functions with different methods for different input types. Below is a chemistry example problem. I do not understand one of the allocations involved.

abstract type RateParams end
 
struct Reaction
    rateparams::RateParams
end

struct elementary_reaction_rate_params <: RateParams
    rate::Array{Float64,1}
end

function reaction_rate!(param::elementary_reaction_rate_params,k::Array{Float64,1})
    for i in 1:10
        k[i] = param.rate[1]*param.rate[2]*param.rate[3]
    end
    nothing
end

rx1 = Reaction(elementary_reaction_rate_params([1.0,2.0,3.0]))

@time k = zeros(Float64,10)
@time rateparams = rx1.rateparams
@time reaction_rate!(rateparams,k)

This returns

  0.000002 seconds (1 allocation: 160 bytes)
  0.000005 seconds (1 allocation: 16 bytes)
  0.000001 seconds

The first allocation for k makes sense. But why is there 1 allocation for the line rateparams = rx1.rateparams? Is there a way to get rid of this allocation? Any other relevant advice is appreciated. Thanks!

[nilshg]

Don’t benchmark in global scope, and use BenchmarkTools.@btime to get accurate benchmarking results:

julia> function test1(rx)
           rateparams = rx.rateparams
       end
test1 (generic function with 1 method)

julia> function test2(rateparams, k)
           reaction_rate!(rateparams, k)
       end
test2 (generic function with 1 method)

julia> @btime test1($rx1)
  1.794 ns (0 allocations: 0 bytes)
elementary_reaction_rate_params([1.0, 2.0, 3.0])

julia> @btime test2($rx1.rateparams, $k)
  24.451 ns (0 allocations: 0 bytes)

[DNF]

@Nicholaswogan

struct Reaction
    rateparams::RateParams
end

For performance reasons you should not use fields with abstract types. You should therefore preferably make this type parametric.

@Nicholaswogan

struct elementary_reaction_rate_params

There is a strong convention that type names should be CapitalCase.

[Nicholaswogan]

Thanks for the advice.

Is it possible to have a function that works quickly on an array of abstract types? For example suppose made the array Array{Reaction,1}, then passed it to the following function

function rates!(reacts::Vector{Reaction},kk::Array{Float64,2})
    for i in 1:length(reacts)
        rx = reacts[i]
        rateparams = rx.rateparams
        reaction_rate!(rateparams,kk[i,:])
    end
    nothing
end

Is it generally bad practice to have arrays of abstract types?

[lmiq]

@Nicholaswogan

Is it generally bad practice to have arrays of abstract types?

Yes, it is. You should avoid using abstract types as much as possible in computations. If each element of reacts there can be of a different type functions will have to be dispatched at run time (selected from a table of methods, for example), and that hurts performance a lot.

Not sure if it was clear from previous answers, but you can do something like this in your example:

julia> struct Reaction{T<:RateParams}
           rateparams::T
       end

julia> struct elementary_reaction_rate_params <: RateParams
           rate::Array{Float64,1}
       end

julia> Reaction(elementary_reaction_rate_params([1.0,2.0,3.0]))
Reaction{elementary_reaction_rate_params}(elementary_reaction_rate_params([1.0, 2.0, 3.0]))

in which case an instance of Reaction is parametrized to the type of rateparams and functions will specialize for this parametric type of Reaction.

You can still have performance problems if you have an array of Reaction instances and each instance is of a different type (if you had, for example, a non_elementary... type, also a subtype of RateParams). The performance problem there will be important if computations done with each element of the array of Reactions are fast in such a way that the delay associated to choosing which method should be used for each element is relevant relative to the time taken on the computations for each element.

[DNF]

@Nicholaswogan

Is it generally bad practice to have arrays of abstract types?

It’s not great for performance, and you should avoid it if you can, but sometimes that’s just not possible, and it’ll have to do.

But you should make sure that you’re not doing it needlessly. And in this case the issue isn’t that you have an array of abstract types, but that you have a concrete type with an abstract field. So you should turn it into a parameterized type, like @lmiq says.

[lmiq]

I profit to make thinks clearer for myself.

Look at this example, which I think illustrates your original problem:

julia> abstract type MyAbstractType end

julia> struct MyType1 <: MyAbstractType
       end

julia> struct MyType2 <: MyAbstractType
       end

julia> mutable struct MyContainerType1
         x::MyAbstractType
       end

I created two types, subtypes of MyAbstractType, and a mutable structure that can contain one element of those types. I can have an instance of MyContainerType1 with:

julia> a = MyContainerType1(MyType1())
MyContainerType1(MyType1())

julia> a.x
MyType1()

and the x element can be changed to MyType2 if I want to:

julia> a.x = MyType2()
MyType2()

julia> a.x
MyType2()

This is flexible, but no method that works on MyContainerType1 will be able to specialize its computations to either of the MyType structs.

On the other side, if I use a parametric container, as:

julia> mutable struct MyContainerType2{T<:MyAbstractType}
         x::T
       end

julia> a = MyContainerType2(MyType1())
MyContainerType2{MyType1}(MyType1())

Note that the type contained is written in the signature of a. When a is passed to a function, that function knows exactly what it is dealing with and can specialize accordingly.

Of course, this comes with a flexibility penalty:

julia> a.x = MyType2()
ERROR: MethodError: Cannot `convert` an object of type MyType2 to an object of type MyType1

You cannot change anymore the type of a.x (even if the structure is mutable, which was intended here). If the structure is immutable you cannot change it either way, but still the dispatch system cannot know in advance if the container type has one or other MyType type inside, and thus cannot specialize either.

[Nicholaswogan]

Thanks so much for the more in depth example. I now understand the importance of using parametric containers.

It seems like, for calculating rates of many reactions, relying on multiple dispatch to switch between reaction types is probably not the most efficient thing because you must operate on a Vector{Reaction}, and the calculation of the rate is very fast. Instead, it is probably better to use if statements without dispatching over a function. Below, I show what I mean.

This makes me a bit sad, because then using ifs reduces the generality of the code.

Using multiple dispatch:

abstract type RateParams end

struct ElementaryReactionRateParams{T<:AbstractFloat} <: RateParams
    a::T
    b::T
    c::T
end

struct ThreeBodyReactionRateParams{T<:AbstractFloat} <: RateParams
    a::T
    b::T
    c::T
end

struct Reaction{T<:RateParams}
    rateparam::T
end

struct ReactionsDispatch{T<:Vector{<:Reaction}}
    rxs::T
end

function reaction_rate(rxparams::ElementaryReactionRateParams)
    return rxparams.a + rxparams.b + rxparams.c
end

function reaction_rate(rxparams::ThreeBodyReactionRateParams)
    return rxparams.a * rxparams.b * rxparams.c
end

function rates_dispatch!(rxs,k::Array{<:AbstractFloat})
    for i in 1:length(rxs.rxs)
        k[i] = reaction_rate(rxs.rxs[i].rateparam)
    end
    nothing
end

# make reaction
rx1 = Reaction(ElementaryReactionRateParams(2.0,2.0,4.0))
rx2 = Reaction(ThreeBodyReactionRateParams(2.0,2.0,4.0))

# now make reactions
rxs = ReactionsDispatch([rx1, rx2]);

using BenchmarkTools
k = Vector{Float64}(undef,length(rxs.rxs))
@btime rates_dispatch!(rxs,k)

result: 170.056 ns (4 allocations: 96 bytes)

Using if statements:

struct ReactionsIf{T<:Array{<:AbstractFloat,2},D<:Vector{<:AbstractString}}
    rxs::T
    rxstype::D
end

# make reaction
rx1 = [2.0 2.0 4.0]
rx2 = [2.0 2.0 4.0]
rx1type = "Eelementary"
rx2type = "ThreeBody"

rxs = [rx1; rx2]
rxstype = [rx1type; rx2type];

reacs = ReactionsIf(rxs,rxstype)

function rates_if!(reacs,k)
    for i in 1:size(reacs.rxs,1)
        if reacs.rxstype[i] == "Eelementary"
            k[i] = reacs.rxs[i,1]*reacs.rxs[i,2]*reacs.rxs[i,3]
        elseif reacs.rxstype[i] == "ThreeBody"
            k[i] = reacs.rxs[i,1]*reacs.rxs[i,2]*reacs.rxs[i,3]
        end
    end
end

using BenchmarkTools
k = Vector{Float64}(undef,size(reacs.rxs,1))
@btime rates_if!(reacs,k)

result: 50.803 ns (0 allocations: 0 bytes)

[lmiq]

Yes, there are quite a few threads here dealing with that type of problem. There is no magic, but one possible solution was given here:

Macro to write function with many conditionals

The minimum working case would be: Let say I have a vector of mixed types: x = [ isodd(i) ? rand(Int64) : rand(Float64) for i in 1:1000 ] I want to sum these numbers, with: function mysum(x) s = 0. for val in x s += val end s end This turns out to the be slow because of the type instability of the numbers and the runtime dispatch associated with it. This version which splits the calculations for each type runs much faster and does not allocate: function mysum_fast(x) s = 0. …