Some doubts about types

[baptiste]

I’ve ported a bit of code from R/Matlab/C++ to Julia a couple of times over the past few years and can get things to work, but I’m still not familiar with the best-practice, particularly around type declarations. My goal is improved efficiency – I already have working code in the 3 languages mentioned above –, but also learning more about Julia as a side-effect.

Here’s a sample of the code that’s giving me some doubts (if you spot other things that seem odd please let me know),

using Base.Threads
using LinearAlgebra

struct Material
    wavelength::Array{Float64}
    alpha::Array{Complex}
    medium::Real
end

function lorentzian(λ, α_k, λ_k, µ_k)
    - α_k*λ_k/µ_k * (1.0 - 1.0 / ( 1.0 - (λ_k/ λ)^2 - 1im*(λ_k^2 /(µ_k*λ))))
end

function alpha_bare(λ::Float64, α_inf::Float64=9.6e-39,
    α_k::Float64=5.76e-38, λ_k::Float64=526.0, µ_k::Float64=1.0e4)

    ε_0 = 8.854e-12
    nm3 = 1e27

    α_trace = α_inf
    for kk in 1:length(α_k)
        α_trace += lorentzian(λ, α_k[kk], λ_k[kk], µ_k[kk])
    end
    prefact = nm3/(4*pi*ε_0)
    prefact * α_trace
end

function polarisation(E::Array{Complex{Float64}}, AlphaBlocks::Array{Complex})

    N = size(AlphaBlocks,3)
    Ni = size(E,2)
    P = similar(E)
    α_tmp = Array{Complex}(undef, (3,3))

    # loop over N particles
    for ii in 1:N
        α_tmp =  AlphaBlocks[:,:,ii]
        ind = ((ii-1)*3 + 1):(ii*3)
        P[ind,1:Ni] = α_tmp * E[ind,1:Ni];
    end

    return P
end

function extinction(kn, P::Array{Complex{Float64}}, Ein::Array{Complex})

    Ni = size(P,2)
    cext = Vector{Float64}(undef, Ni)

    for jj in 1:Ni
        cext[jj] = imag(dot(Ein[:,jj], P[:,jj])) # E^* P
    end

    return  4*pi*kn*cext
end

function euler_active(φ::Float64, θ::Float64, ψ::Float64)

    R = Array{Float64}(undef,(3,3))
    cosφ = cos(φ); cosψ = cos(ψ); cosθ = cos(θ)
    sinφ = sin(φ); sinψ = sin(ψ); sinθ = sin(θ)

    R[1,1] = cosφ*cosθ*cosψ - sinφ*sinψ
    R[2,1] = sinφ*cosθ*cosψ + cosφ*sinψ
    R[3,1] = -sinθ*cosψ

    R[1,2] = -cosφ*cosθ*sinψ - sinφ*cosψ
    R[2,2] = -sinφ*cosθ*sinψ + cosφ*cosψ
    R[3,2] = sinθ*sinψ

    R[1,3] = cosφ*sinθ
    R[2,3] = sinφ*sinθ
    R[3,3] = cosθ

    return R

end

The core calculation deals with linear algebra of real and complex matrices, and I define a couple of custom types to wrap such objects in a convenient structure.

Some questions about the above: I’ve used things like Array{Float64}(undef,(3,3)) to initialise arrays, and P::Array{Complex{Float64}} to specify complex types in function arguments. Are these the recommended types here? I thought I should switch to the most general type where the code might make sense, such as replacing Float64 with Real throughout, but my code stopped working and I wasn’t able to identify or even locate the error. I find it all the more difficult to diagnose the problem that my custom structs cannot be redefined within a session (or can they?), so every time I change the types inside I need to restart.

Thanks for any pointers.

[stevengj]

I think you probably want something like:

struct Material{T<:Real}
    wavelength::Vector{T}
    alpha::Vector{Complex{T}}
    medium::T
end

That way, you aren’t defining a single type, you are defining a family of types Material{T} for any real type that you might want, e.g. Material{Float64} for double precision. This gives allows you to write both generic code (works for any T) while still retaining the performance advantage of concrete types. (In contrast, Array{Complex} is an array of abstract types, i.e. an array of pointers to “boxed” objects with type tags that will be unbearably slow, and Array is also abstract because it can be any dimensionality. Vector{T} is an alias for the concrete type Array{T,1}, i.e. a 1d array of T.)

For example, see the tutorial here on the utility and performance of user-defined parameterized types.

You could go even further, and allow a wider variety of types for your fields, e.g. different array types and different real types for the different fields, by adding more type parameters:

struct Material{TW<:AbstractVector{<:Real}, TA<:AbstractVector{<:Complex}, TM<:Real}
    wavelength::TW
    alpha::TA
    medium::TM
end

However, putting in more genericity than you need can lead to painfully parameterized types, so it’s a tradeoff in coding (but not in performance).

[baptiste]

That’s very helpful, thanks!

The different kinds of arrays are stlll a bit mysterious to me, between Vector, Matrix, Array, AbstractArray, RowVector(which I saw mentioned but doesn’t seem to exist).
I’m inferring that Vector is an alias for 1D array, and probably Matrix an alias for 2D arrays. This leaves me with the AbstractArray concept – probably overly general for a regular kind of array?

And then to initialise or define types of their entries I’m not sure if I should go for Complex{Float64}, or Complex{AbstractFloat}, or Complex{Real}.

The {} notation remains a bit of a mystery for me, especially when nested, as in Array{Complex{Float64}}(undef, (3,3)); I wonder if there’s a dedicated learning resource for these, focusing on linear algebra a la Matlab?

Below’s an additional excerpt from the code where my type definitions may also be suboptimal,

function propagator_freespace_labframe(R::Matrix{Float64}, kn::Float64, AlphaBlocks::Array{Complex{Float64}})
    N = size(R,2)
    A = Matrix{Complex{Float64}}(I, 3*N, 3*N)
    I3 = Matrix{Complex{Float64}}(I, 3, 3)

    rk_to_rj = Vector{AbstractFloat}(undef,3)
    rjkhat = Vector{AbstractFloat}(undef,3)
    rjkrjk = Array{AbstractFloat}(undef,(3, 3))
    Ajk = Array{Complex{Float64}}(undef,(3, 3))
    α_jj = similar(Ajk)
    α_kk = similar(Ajk)

    # nested for loop over N dipoles
    for jj in 1:N

        α_jj =  AlphaBlocks[:,:,jj]

        for kk in (jj+1):N
            α_kk =  AlphaBlocks[:,:,kk]

            rk_to_rj = R[:,jj] - R[:,kk]
            rjk = norm(rk_to_rj, 2)
            rjkhat = rk_to_rj / rjk
            rjkrjk =  rjkhat * transpose(rjkhat)

            Ajk = exp(im*kn*rjk) / rjk *  (kn*kn*(rjkrjk - I3) +
            (im*kn*rjk - 1.0) / (rjk*rjk) * (3*rjkrjk - I3)) ;

            # assign block
            ind1 = ((jj-1)*3 + 1):(jj*3)
            ind2 = ((kk-1)*3 + 1):(kk*3)
            A[ind1, ind2] = Ajk * α_kk
            # symmetric block
            A[ind2, ind1] = transpose(Ajk) * α_jj

        end
    end

    return A
end

Edit: I’ve run into issues with Matrix when the number of columns happens to be just 1; I guess a Vector is a subtype of Array, but not of a Matrix? I’ll stick with Array I think.

(PS: you’ll recognise obviously a coupled-dipole implementation)

[kristoffer.carlsson]

@baptiste

And then to initialise or define types of their entries I’m not sure if I should go for Complex{Float64}, or Complex{AbstractFloat}, or Complex{Real}.

Here, go with Complex{Float64}. But instead of initializing as

Ajk = Array{Complex{Float64}}(undef,(3, 3))

it is much easier to use something like fill(0.0 + 0.0im, 5, 5) or zeros(ComplexF64, 5, 5). Note that ComplexF64 is a short way of writing Complex{Float64}.

You would want to use Array{Complex{Real}} if you want to be able to store a bunch of numbers with different types (all being a subtype of Read) in the array, but this is quite rare.

@baptiste

I guess a Vector is a subtype of Array, but not of a Matrix?

Vector is not a subtype of Array, it is an Array with 1 dimension, i.e. Vector{T} == Array{T, 1}. And Matrix{T} == Array{T, 2}.

[stevengj]

@baptiste

I’m inferring that Vector is an alias for 1D array, and probably Matrix an alias for 2D arrays. This leaves me with the AbstractArray concept – probably overly general for a regular kind of array?

Array is the built-in multidimensional array type, and in particular it is a family of types Array{T,N}: N-dimensional arrays of element type T. Vector{T} is an alias for Array{T,1} and Matrix{T} is an alias for Array{T,2}.

Array is a subtype of the abstract AbstractArray type — by accepting AbstractArray rather than Array, you allow your code to be used with other array types. An example in Base is a range like 1:100, which is acts like a 1d array but its elements are computed as needed. Other examples in external packages include StaticArrays and GPUArrays.

In general, the more people you expect to re-use your code, the more generic you should try to be, so that they can combine your code arbitrarily with different packages. But as long as it is mostly you using your code, you can always go back and add more type parameters to make your code more generic. As you get more experience working with Julia, you should learn to instinctively write more generic code by default.

@baptiste

And then to initialise or define types of their entries I’m not sure if I should go for Complex{Float64}, or Complex{AbstractFloat}, or Complex{Real}.

You almost always want to parameterize by concrete types. For an explanation of why this is a critical distinction, see:

Problem with Complex{Rationals}

You should also understand that this distinction is absolutely critical for performance. It’s not just an arbitrary pedantic choice. Think about how instances of the types are actually stored and processed on the computer. A Complex{Rational{Int}} for a+bi can be (and is) stored as four consecutive Int (64-bit integer) values: (anum,aden,bnum,bden). For operations on such quantities, since the compiler knows everything’s type and how it is stored, it can load the numerators and denominators …

[stevengj]

@baptiste

function propagator_freespace_labframe(R::Matrix{Float64}, kn::Float64, AlphaBlocks::Array{Complex{Float64}}) N = size(R,2) A = Matrix{Complex{Float64}}(I, 3*N, 3*N) I3 = Matrix{Complex{Float64}}(I, 3, 3)

For function arguments, you almost always want to be more generic. Indeed, you can leave off the type declarations completely and it is perfectly fine — Julia will type-specialize your function at compile time based on the type of arguments that you pass. But it can be nice to add type declarations to indicate the expected types, for clarity, and also to allow you to define different methods for different types.

For example, a very generic definition of your function would be something like:

function propagator_freespace_labframe(R::AbstractMatrix{<:Real}, kn::Real, AlphaBlocks::AbstractArray{<:Complex,3})
    N = size(R,2)
    T = float(promote_type(eltype(R), typeof(kn), real(eltype(AlphaBlocks))))
    A = Matrix{Complex{T}}(I, 3*N, 3*N)
    I3 = Matrix{Complex{T}}(I, 3, 3)

    rk_to_rj = Vector{T}(undef,3)
    rjkhat = Vector{T}(undef,3)
    rjkrjk = Array{T}(undef,(3, 3))
    Ajk = Array{Complex{T}}(undef,(3, 3))

This way, you accept any type of argument, and the type of the temporary arrays uses a floating-point type T that is chosen based on the types that are passed as arguments.

On a separate note, note that you are allocating tons of copies of your arrays with slicing expressions like AlphaBlocks[:,:,jj], and your vectorized expressions allocate zillions of temporary arrays. You need to be cautious about writing Matlab-style code in Julia, which does not take full advantage of the language. See also the More Dots: Syntactic Loop Fusion in Julia blog post.

In this case, it looks like your vectorized operations in your inner loops are operating on 3-component coordinate vectors. For this application, it will be much faster and more elegant to use StaticArrays, in which the 3 dimensions are known at compile time and all of the 3-element loops are inlined.

That is, I’m guessing that AlphaBlocks should be a Vector{SMatrix{3,3,Complex{T}}}, R should be a Vector{SVector{3,T}}, A should be a Matrix{SMatrix{3,3,Complex{T}}}, and so on. Then you would do things like rk_to_rj = R[jj] - R[kk] and A[jj,kk] = Ajk * AlphaBlocks[kk].

Note also that explicitly forming I3 is a waste — the whole point of the I object is that things like rjkrjk - I are computed efficiently by only subtracting the diagonal elements, with I implicitly “expanding” to a matrix of the correct size.

[baptiste]

Thanks so much for all the helpful directions.
The code is indeed a direct port from Matlab (and R, and C++ Armadillo, all 3 with the same basic syntax), and while I was pretty sure this wouldn’t produce good Julia code, it’s hard for a novice to know the “rights” and “wrongs” without some sort of code review by someone experienced in the language.
I’ll rewrite it with these new directions and hopefully will get closer to internalising the language’s philosophy with each iteration.

[baptiste]

@stevengj

On a separate note, note that you are allocating tons of copies of your arrays with slicing expressions like AlphaBlocks[:,:,jj] , and your vectorized expressions allocate zillions of temporary arrays. You need to be cautious about writing Matlab-style code in Julia, which does not take full advantage of the language. See also the More Dots: Syntactic Loop Fusion in Julia blog post.

Thanks for the hints, it’s definitely a case of forcing oneself to unlearn Matlab/R patterns; I’ve tried to rewrite that piece of code following your tips and it now looks like this,

function propagator_freespace_labframe(R::Vector{SVector{3,Float64}},
    kn::Float64,
    Alpha::Vector{SVector{3,Complex{Float64}}},
    Angles::Vector{SVector{3,Float64}})

    N = length(R)
    
    A = Matrix{SMatrix{3,3,Complex{Float64}}}(undef, N, N)
    eye3 = SMatrix{3,3}(Array{Complex{Float64}}(I, 3, 3))
    Ajk = MMatrix{3,3,Complex{Float64}}
    α_jj = similar(Ajk)
    α_kk = similar(Ajk)

    # nested for loop over N dipoles
    for jj in 1:N

        A[jj,jj] = eye3

        Rm = euler_passive(Angles[jj]...)
        α_jj = transpose(Rm) * Diagonal(Alpha[jj]) * Rm


        for kk in (jj+1):N

            Rm = euler_passive(Angles[kk]...)
            α_kk = transpose(Rm) * Diagonal(Alpha[kk]) * Rm

            rk_to_rj = R[jj] - R[kk]

            rjk = norm(rk_to_rj, 2)
            rjkhat = rk_to_rj / rjk
            rjkrjk =  rjkhat * transpose(rjkhat)

            Ajk = SMatrix(exp(im*kn*rjk) / rjk *  (kn*kn*(rjkrjk - I) +
            (im*kn*rjk - 1.0) / (rjk*rjk) * (3*rjkrjk - I)))

            # assign block
            A[jj,kk] = Ajk * α_kk
            # symmetric block
            A[kk,jj] = transpose(Ajk) * α_jj

        end
    end

    return A
end

# testing with dummy values
R = [@SVector randn(3) for ii=1:4]
Angles = [@SVector randn(3) for ii=1:4]
kn = 0.2
Alpha = [@SVector randn(Complex{Float64},3) for ii=1:4]

propagator_freespace_labframe(R,kn,Alpha,Angles)

Is this more or less along the lines of what you suggested? I still explicitly expanded I for the diagonal blocks of A, otherwise with just I I got an error, but in Ajk I now use I instead of the expanded 3x3 identity.
My biggest doubt about this approach though, is that ultimately A appears in a 3Nx3N linear system A x = b (coupled-dipole equations), and if A is a NxN matrix of SMatrices I don’t really see how it’s going to work (ultimately some Lapack routine will need the 3Nx3N matrix, won’t it?). Does it make sense to build A this way, as N^2 blocks (which is definitely more readable), and then somehow convert it to a 3Nx3N matrix?

I’ve read the “more dots” post many times (in fact it’s pushed me to learn more about julia, so thanks!) but I’m not sure how it would apply in this particular function, where I’m building 3x3 blocks. I can’t see a way to combine those lines into one expression,

            rk_to_rj = R[jj] - R[kk]

            rjk = norm(rk_to_rj, 2)
            rjkhat = rk_to_rj / rjk
            rjkrjk =  rjkhat * transpose(rjkhat)

            Ajk = SMatrix(exp(im*kn*rjk) / rjk *  (kn*kn*(rjkrjk - I) +
            (im*kn*rjk - 1.0) / (rjk*rjk) * (3*rjkrjk - I)))

[stevengj]

@baptiste

My biggest doubt about this approach though, is that ultimately A appears in a 3Nx3N linear system A x = b (coupled-dipole equations), and if A is a NxN matrix of SMatrices I don’t really see how it’s going to work (ultimately some Lapack routine will need the 3Nx3N matrix, won’t it?).

You can do Ajk = ...some SMatrix.. and then do A[3j-2:3j, 3k-3:3k] = Ajk to write the Smatrix into the “flat” array A in-place.

@baptiste

α_jj = transpose(Rm) * Diagonal(Alpha[jj]) * Rm

You can do this by dot(Rm, Alpha[jj] .* Rm) or Rm' * (Alpha[jj] .* Rm)without constructing the Diagonal matrix.

Note also that you can just do foo' instead of transpose(foo) — the two are equivalent for real vectors like your coordinate vectors here.

@baptiste

Ajk = SMatrix(...)

If you’ve done things correctly, you shouldn’t need the SMatrix(...) constructor here… the matrix expression should already give an SMatrix (because it’s built out of SMatrix expressions like rjkhat * rjkhat').

(You can add a type assertion ::SMatrix at the end of the expression if you want to check that you did this correctly. A type assertion will throw an error if the expression is not the expected type, but is eliminated at compile-time if the expression has the correct type and type-inference succeeds.)