How to specialize arguments to avoid unions

New to Julia
,
1

Hello everyone!

I am trying to write two functions that accept 3D and 2D arrays that can be either real or complex. The first function receives arguments of one dimension higher than the second function. Then, it iterates over that extra dimension calling the second function for each “slice” of the arguments along that dimension. The purpose of the code is converting S-matrices for many frequencies at given reference impedances for all ports, also at that many frequencies.

The problem is that all arguments can be either Real or Complex arrays and even different between them, although the most common case is that the first argument is complex and the rest are real. How can I write the header of the function to admit Real and Complex arrays without explicitly typing

CR = Union{Real,Complex}
function S2S(S::Array{T1,3},Za::Array{T2,2},Zb::Array{T3,2}) where {T1 <:CR && T2 <: CR && T3 <: CR}
...
end

I leave my current code here. Thanks!

function S2S(Sₐ::Array{ComplexF64,3},Zₐ::Array{Float64,2},Zᵦ::Array{Float64,2})
    Sᵦ = similar(Sₐ)
    for f in 1:size(Sₐ,3)
        Sᵦ[:,:,f] = S2S(Sₐ[:,:,f],Zₐ[:,f],Zᵦ[:,f])
    end
    return Sᵦ
end

function S2S(Sₐ::Array{ComplexF64,2},Zₐ::Array{Float64,1},Zᵦ::Array{Float64,1})
    U = UniformScaling(size(Sₐ,1))
    G = diagm(sqrt.(Zₐ./Zᵦ))
    Γ = diagm((Zᵦ-Zₐ)./(Zᵦ+Zₐ))
    return Sᵦ = G*inv(U-Sₐ)*(Sₐ-Γ)*inv(U-Sₐ*Γ)*(U-Sₐ)*inv(G)
end
2

Do you really need those type annotations? Why not AbstractArray everywhere, or even not annotating at all?

Lyndon White has a nice blog post on “Over-constraining argument types”, it is one of the sections of the post here. Maybe that helps.

3

I do need to type the dimensions so the first function calls the second. I can substitute AbstractArray where Array is typed, but otherwise, I cannot see how this would help.

Why is better AbstractArray than Array?

Thanks!

4

It will accept @views, for example.

Edit: something like this will be more concise, anyway:

julia> function g!(x::AbstractArray{T1,3},y::AbstractArray{T2,2}) where {T1,T2 <: Number}
         x = 2*x
         y = 2*y
         nothing
       end
g! (generic function with 1 method)

julia> x = rand(ComplexF64,3,3,3);

julia> y = rand(3,2);

julia> g!(x,y)

(adapted to your case, of course)

2 Likes
5

With regards to choosing AbstractArray rather than Array, consider that there are many subtypes of AbstractArray that may not appear relevant when first writing your algorithm, but that generalizing to is very convenient (Adjoint, UniformScaling, Diagonal, StaticArray just to name a few!).

If the first argument must be a complex number (that is, real arguments are completely not allowed ever), then it makes sense to restrict it. Otherwise, go with Number (which already is the supertype of Real and Complex, so no need to define CR). For the other arguments, restricting to Float64 is almost never a good idea in practice. Should your algorithm really fail if the inputs are integers, or would all of the math “just work” anyway? I would probably consider a signature like this (again, assuming complex is absolutely required for the first one, otherwise I would leave that eltype off as well)

function S2S(Sₐ::AbstractArray{<:Complex, 3}, Zₐ::AbstractMatrix, Zᵦ::AbstractMatrix)
    ...
end

function S2S(Sₐ::AbstractMatrix{<:Complex}, Zₐ::AbstractVector, Zᵦ::AbstractVector)
    ...
end
7 Likes
6

Thanks, super helpful!!

Just the last question: why is it better to have AbstractMatrix than AbstractMatrix{Number}? Shouldn’t it be better this other way, in order to prevent, for example, matrices of String?

7

Parametric types are not covariant. Using Matrix{Number} will not match an argument of Matrix{Float64}, for instance. Instead, you should specify Matrix{<:Number}.

8

Yeah, sorry, I meant to put the <: but forgot it.

In any case, is it better for the compiler not to specify the least possible type, or there is no usage and performance penalty?

9

The annotations of types on function arguments are only for dispatch. It makes no difference whatsoever for performance. The function will be specialized to the given type when executed. This confusion is one of the reasons that make we tend to start annotating everything (the blog post of Lyndon White above is very good in emphasizing that - it discusses the other reasons for using annotations, of course, but performance is not one of those).

2 Likes
10

The usage penalty is that you may constrain it to be too narrow and disallow some valid use cases.

4 Likes