Performance differences when looping over arrays of structs

TLDR version: Vector{Gaussian} should be Vector{Gaussian{Float32}}.

Your struct definition Gaussian{T<:Float32} is somewhat nonsensical. Float32 cannot have any subtypes, so T can only be a Float32. Therefore, either remove the restriction, replace it with something else (e.g. Real), or just take away the type parameter entirely.

This is bad:

gaussians::Array{Gaussian} = []

It creates a Vector{Any} which it then converts to a Vector{Gaussian}. Worse, Vector{Gaussian} is not a Vector{Gaussian{Float32}}. So either remove the whole type parameter, or make sure to use it. So,

# good:
gaussians = Vector{Gaussian{Float32}}()
gaussians = Gaussian{Float32}[] # same as above

# bad
gaussians = Vector{Gaussian}()
# very bad, don't use this style, put types on the right hand side when constructing.
gaussians::Vector{Gaussian} = []

Same here, bad style

total::Float32 = 0.

Do this instead

total = Float32(0.0)
# or use Float32 literal
total = 0.0f0
# or the generic way
total = zero(Float32) 

In dynamic languages, types belong to values, not to variables.

BTW, you'll have to modify some of your function definitions:

total_intensity1(gaussian_list::Array{Gaussian})

should be

total_intensity1(gaussian_list::Array{<:Gaussian})

There's more, but this is a start.

Edit: OK, a few more things:

1. rand(Float32, (length, 1)) length is a super important function in Base, so it's normally good not to shadow it like this. And, make vectors instead of matrices:

   rand(Float32, (N, 1)) # this is an Nx1 matrix

   rand(Float32, N) # this is a length-N vector

2. push!(gaussians, new_gaussian) This iteratively resizes the vector over and over. When you know the size of the vector as in your case, it is better to pre-allocate:

   gaussians = Vector{Gaussian{Float32}}(undef, 2000)

3. You can avoid an unnecessary allocation here:

   total = sum(intensity.(gaussian_list, Float32(0.75), Float32(0.11))) like this:

   total = sum(g->intensity(g, 0.75f0, 0.11f0), gaussian_list)

Explanation: sum(f.(x)) first creates the array f.(x), then sums it, while sum(f, x) just applies f to each element before adding it to the sum.

Here's an implementation with benchmarks:

struct Gaussian{T<:Real}
    x0::T
    y0::T
    A::T
    a::T
    b::T
    c::T
end
Gaussian(x::Real...) = Gaussian(promote(x...)...)

function intensity(model::Gaussian, x::Real, y::Real)
    val = model.A * exp(
        -( 
            model.a * (x - model.x0)^2 +
            2 * model.b * (x - model.x0) * (y - model.y0) +
            model.c * (y - model.y0)^2
        )
    )
    return val
end

function build_array(N=2000)
    return [Gaussian(ntuple(_->rand(Float32), 6)...) for _ in 1:N]
end

Benchmarks (remember to interpolate variables, and avoid global scope):

julia> gaussians = build_array(2000);

julia> @btime sum(intensity.($gaussians, Float32(0.75), Float32(0.11)))
  14.600 μs (1 allocation: 7.94 KiB)
947.5305f0

julia> @btime sum(g->intensity(g, 0.75f0, 0.11f0), $gaussians)
  12.600 μs (0 allocations: 0 bytes)
947.5309f0

There's a slight difference in the final sums, since sum of vector uses a numerically superior method called pairwise summation.

Final bonus: Try Tullio.jl, which also uses multithreading. It doesn't make any difference for 2000 elements, but it does for larger arrays (using 12 threads here):

julia> using Tullio, LoopVectorization

julia> gaussians = build_array(200_000);

julia> @btime sum(g->intensity(g, 0.75f0, 0.11f0), $gaussians)
  1.228 ms (0 allocations: 0 bytes)
92722.7f0

julia> @btime @tullio s := intensity($gaussians[i], 0.75f0, 0.11f0)
  330.100 μs (197 allocations: 11.53 KiB)
92722.79f0

I don't have time to figure this one out in detail unfortunately, but the first thing I'd say is check whether this is a benchmarking artifact - gaussians is a global variable which should be interpolated into the benchmark using $.

As to your function, the type annotations are not doing anything for performance here, and will make your function less composable (e.g. you won't be able to autodiff through it give you're restricting everything to Float32).

Here's how I would write it:

function intensity(m, x, y)
    (; x₀, y₀, A, a, b, c) = m # destructuring input
    A * exp( -(a * (x - x₀)^2 + 2b * (x - x₀) * (y - y₀) + c * (y - y₀)^2 ) )
end

With that I'm getting:

  231.100 μs (12001 allocations: 195.44 KiB)
  231.500 μs (12001 allocations: 195.44 KiB)
  229.200 μs (12000 allocations: 187.50 KiB)
  229.300 μs (12000 allocations: 187.50 KiB)

which is about 100μs faster than your original version on my machine.