# Some doubts about types

**URL:** <https://discourse.julialang.org/t/some-doubts-about-types/27635>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [August 16, 2019, 11:05pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635 "2019-08-16T23:05:41Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![baptiste](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptiste/32/9813_2.png) [@baptiste](https://discourse.julialang.org/u/baptiste)\
**Post date:** [August 16, 2019, 11:05pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/1 "2019-08-16T23:05:41Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 17, 2019, 2:38am UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/2 "2019-08-17T02:38:15Z")

</div>

I think you probably want something like:

```julia
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](https://nbviewer.jupyter.org/github/mitmath/18S096/blob/master/lectures/lecture2/Code-specialization.ipynb#Defining-our-own-types) 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:

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

---

<div class="post-metadata">

**Author:** ![baptiste](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptiste/32/9813_2.png) [@baptiste](https://discourse.julialang.org/u/baptiste)\
**Post date:** [August 17, 2019, 3:38am UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/3 "2019-08-17T03:38:14Z")

</div>

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,

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

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [August 17, 2019, 7:18am UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/4 "2019-08-17T07:18:31Z")

</div>

> [@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

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 17, 2019, 1:18pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/5 "2019-08-17T13:18:32Z")

</div>

> [@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}](https://discourse.julialang.org/t/problem-with-complex-rationals/9474/7):
>
> 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 …

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 17, 2019, 1:39pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/6 "2019-08-17T13:39:26Z")

</div>

> [@baptiste](#):
>
> ```julia
> 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:

```julia
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](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](https://github.com/JuliaArrays/StaticArrays.jl), 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.

---

<div class="post-metadata">

**Author:** ![baptiste](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptiste/32/9813_2.png) [@baptiste](https://discourse.julialang.org/u/baptiste)\
**Post date:** [August 17, 2019, 9:39pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/7 "2019-08-17T21:39:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 18, 2019, 2:38pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/8 "2019-08-18T14:38:55Z")

</div>

7 posts were split to a new topic: [Grassman algebras for electromagnetic scattering](https://discourse.julialang.org/t/grassman-algebras-for-electromagnetic-scattering/27681)

---

<div class="post-metadata">

**Author:** ![baptiste](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baptiste/32/9813_2.png) [@baptiste](https://discourse.julialang.org/u/baptiste)\
**Post date:** [August 31, 2019, 9:40am UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/9 "2019-08-31T09:40:33Z")

</div>

> [@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](https://julialang.org/blog/2017/01/moredots) 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,

```julia

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,

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

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [August 31, 2019, 2:00pm UTC](https://discourse.julialang.org/t/some-doubts-about-types/27635/10 "2019-08-31T14:00:00Z")

</div>

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