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 `struct`s 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.

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).

8 Likes

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)

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.

`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}`.

6 Likes

`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.

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

8 Likes

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 https://julialang.org/blog/2017/01/moredots 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.

8 Likes

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.

7 posts were split to a new topic: Grassman algebras for electromagnetic scattering

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)))
``````

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.

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.

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.)

2 Likes