Skip to content

The Importance of Type-Stability

You are here because you have probably heard about the excellent performance of Julia compared to other common programming languages like Python. One of the reasons is the Just-In-Time (JIT) compiler of Julia, which is able to generate highly optimized machine code. However, the JIT compiler can only do its job if the code type can be inferred. You can also read the Performance Tips section in Julia's documentation for more details. Here, we try to explain it briefly, with a focus on the QuantumToolbox.jl package.

Note

This page is not a tutorial on QuantumToolbox.jl, but rather a general guide to writing Julia code for simulating quantum systems efficiently. If you don't care about the performance of your code, you can skip this page.

Basics of type stability

Let's have a look at the following example:

julia
function foo(x)
    if x > 0
        return 1
    else
        return -1.0
    end
end

The function foo apparently seems to be innocent. It takes an argument x and returns either 1 or -1.0 depending on the sign of x. However, the return type of foo is not clear. If x is positive, the return type is Int, otherwise it is Float64. This is a problem for the JIT compiler, because it has to determine the return type of foo at runtime. This is called type instability (even though it is a weak form) and may lead to a significant performance penalty. To avoid this, always aim for type-stable code. This means that the return type of a function should be clear from the types of its arguments. We can check the inferred return type of foo using the @code_warntype macro:

julia
@code_warntype foo(1)
MethodInstance for Main.var"Main".foo(::Int64)
  from foo(x) @ Main.var"Main" type_stability.md:18
Arguments
  #self#::Core.Const(Main.var"Main".foo)
  x::Int64
Body::UNION{FLOAT64, INT64}
1 ─ %1 = Main.var"Main".:>::Core.Const(>)
│   %2 = (%1)(x, 0)::Bool
└──      goto #3 if not %2
2 ─      return 1
3 ─      return -1.0

The key point is to ensure the return type of a function is clear from the types of its arguments. There are several ways to achieve this, and the best approach depends on the specific problem. For example, one can use the same return type:

julia
function foo(x)
    if x > 0
        return 1.0
    else
        return -1.0
    end
end

Or you can ensure the return type matches the type of the argument:

julia
function foo(x::T) where T
    if x > 0
        return T(1)
    else
        return -T(1)
    end
end

The latter example is very important because it takes advantage of Julia's multiple dispatch, which is one of the most powerful features of the language. Depending on the type T of the argument x, the Julia compiler generates a specialized version of foo that is optimized for this type. If the input type is an Int64, the return type is Int64, if x is a Float64, the return type is Float64, and so on.

julia
@show foo(1)
@show foo(-4.4)
@show foo(1//2)
foo(1) = 1
foo(-4.4) = -1.0
foo(1 // 2) = 1//1

Note

If you didn't know how to make this function type-stable, it is probably a good idea to read the official Julia documentation, and in particular its Performance Tips section.

Global variables

Another source of type instability is the use of global variables. In general, it is a good idea to declare global variables as const to ensure their type is fixed for the entire program. For example, consider the following function that internally takes a global variable y:

julia
y = 2.4

function bar(x)
    res = zero(x) # this returns the zero of the same type of x
    for i in 1:1000
        res += y * x
    end
    return res
end

The Julia compiler cannot infer the type of res because it depends on the type of y, which is a global variable that can change at any time of the program. We can check it using the @code_warntype macro:

julia
@code_warntype bar(3.2)
MethodInstance for Main.var"Main".bar(::Float64)
  from bar(x) @ Main.var"Main" type_stability.md:79
Arguments
  #self#::Core.Const(Main.var"Main".bar)
  x::Float64
Locals
  @_3::UNION{NOTHING, TUPLE{INT64, INT64}}
  res::ANY
  i::Int64
Body::ANY
1 ─ %1  = Main.var"Main".zero::Core.Const(zero)
│         (res = (%1)(x))
│   %3  = Main.var"Main".:(:)::Core.Const(Colon())
│   %4  = (%3)(1, 1000)::Core.Const(1:1000)
│         (@_3 = Base.iterate(%4))
│   %6  = @_3::Core.Const((1, 1))
│   %7  = (%6 === nothing)::Core.Const(false)
│   %8  = Base.not_int(%7)::Core.Const(true)
└──       goto #4 if not %8
2 ┄ %10 = @_3::Tuple{Int64, Int64}
│         (i = Core.getfield(%10, 1))
│   %12 = Core.getfield(%10, 2)::Int64
│   %13 = Main.var"Main".:+::Core.Const(+)
│   %14 = res::ANY
│   %15 = Main.var"Main".:*::Core.Const(*)
│   %16 = (%15)(Main.var"Main".y, x)::ANY
│         (res = (%13)(%14, %16))
│         (@_3 = Base.iterate(%4, %12))
│   %19 = @_3::UNION{NOTHING, TUPLE{INT64, INT64}}
│   %20 = (%19 === nothing)::Bool
│   %21 = Base.not_int(%20)::Bool
└──       goto #4 if not %21
3 ─       goto #2
4 ┄ %24 = res::ANY
└──       return %24

While in the last example of the foo function we got a weak form of type instability, returning a Union{Int, Float64}, in this case the return type of bar is Any, meaning that the compiler doesn't know anything about the return type. Thus, this function has nothing different from a dynamically typed language like Python. We can benchmark the performance of bar using the BenchmarkTools.jl package:

julia
using BenchmarkTools

@benchmark bar(3.2)
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  40.525 μs …  31.912 ms  ┊ GC (min … max): 0.00% … 99.73%
 Time  (median):     60.011 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   60.134 μs ± 321.203 μs  ┊ GC (mean ± σ):  5.29% ±  1.00%

      ▂                               ▂▂ █▇                     
  ▃▆▂██▄▂▂▂▂▂▁▂▁▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▄▄▃██▆██▇▃▂▂▂▂▂▁▁▂▂▂▂▃▂▂▃▃▂ ▃
  40.5 μs         Histogram: frequency by time         71.4 μs <

 Memory estimate: 46.88 KiB, allocs estimate: 3000.

Here we see a lot of memory allocations and low performances in general. To fix this, we can declare a const (constant) variable instead:

julia
const z = 2.4

function bar(x)
    res = zero(x) # this returns the zero of the same type of x
    for i in 1:1000
        res += z * x
    end
    return res
end

@benchmark bar(3.2)
BenchmarkTools.Trial: 10000 samples with 57 evaluations per sample.
 Range (min … max):  868.789 ns …  3.824 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     869.684 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   876.998 ns ± 41.257 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █     ▂                                                 ▁▁   ▁
  █▄▄▆▄▅█▇▄▄▁▄▃▅▆▇▆▅▁▄▁▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▇████▇ █
  869 ns        Histogram: log(frequency) by time       999 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And we can see that the performance has improved significantly. Hence, we highly recommend using global variables as const, but only when truly necessary. This choice is problem-dependent, but in the case of QuantumToolbox.jl, this can be applied for example in the case of defining the Hilbert space dimensions, static parameters, or the system operators.

Although it is always a good practice to avoid such kind of type instabilities, in the actual implementation of QuantumToolbox.jl (where we mainly deal with linear algebra operations), the compiler may perform only a few runtime dispatches, and the performance penalty may be negligible compared to the heavy linear algebra operations.

Vectors vs Tuples vs StaticArrays

Julia has many ways to represent arrays or lists of general objects. The most common are Vectors and Tuples. The former is a dynamic array that can change its size at runtime, while the latter is a fixed-size array that is immutable, and where the type of each element is already known at compile time. For example:

julia
v1 = [1, 2, 3] # Vector of Int64
v2 = [1.0 + 2.0im, 3.0 + 4.0im] # Vector of ComplexF64
v3 = [1, "ciao", 3.0] # Vector of Any

t1 = (1, 2, 3) # Tuple of {Int64, Int64, Int64}
t2 = (1.0 + 2.0im, 3.0 + 4.0im) # Tuple of {ComplexF64, ComplexF64}
t3 = (1, "ciao", 3.0) # Tuple of {Int64, String, Float64}

@show typeof(v1)
@show typeof(v2)
@show typeof(v3)
@show typeof(t1)
@show typeof(t2)
@show typeof(t3)
typeof(v1) = Vector{Int64}
typeof(v2) = Vector{ComplexF64}
typeof(v3) = Vector{Any}
typeof(t1) = Tuple{Int64, Int64, Int64}
typeof(t2) = Tuple{ComplexF64, ComplexF64}
typeof(t3) = Tuple{Int64, String, Float64}

Thus, we highly recommend using Vector only when we are sure that it contains elements of the same type, and only when we don't need to know its size at compile time. On the other hand, Tuples are less flexible but more efficient in terms of performance. A third option is to use the SVector type from the StaticArrays.jl package. This is similar to Vector, where the elements should have the same type, but it is fixed-size and immutable. One may ask when it is necessary to know the array size at compile time. A practical example is the case of ptrace, where it internally reshapes the quantum state into a tensor whose dimensions depend on the number of subsystems. We will see this in more detail in the next section.

The QuantumObject internal structure

Before making a practical example, let's see the internal structure of the QuantumObject type. As an example, we consider the case of three qubits, and we study the internal structure of the σ^x(2) operator:

julia
σx_2 = tensor(qeye(2), sigmax(), qeye(2))

Quantum Object:   type=Operator   dims=[2, 2, 2]   size=(8, 8)   ishermitian=true
8×8 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries:
     ⋅          ⋅      1.0+0.0im  …      ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅    
 1.0+0.0im      ⋅          ⋅             ⋅          ⋅          ⋅    
     ⋅      1.0+0.0im      ⋅             ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅             ⋅      1.0+0.0im      ⋅    
     ⋅          ⋅          ⋅      …      ⋅          ⋅      1.0+0.0im
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅         1.0+0.0im      ⋅          ⋅

and its type is

julia
obj_type = typeof(σx_2)
QuantumObject{OperatorQuantumObject, Dimensions{3, Tuple{Space, Space, Space}}, SparseMatrixCSC{ComplexF64, Int64}}

This is exactly what the Julia compiler sees: it is a QuantumObject, composed by a field of type SparseMatrixCSC{ComplexF64, Int64} (i.e., the 8x8 matrix containing the Pauli matrix, tensored with the identity matrices of the other two qubits). Then, we can also see that it is a OperatorQuantumObject, with 3 subsystems in total. Hence, just looking at the type of the object, the compiler has all the information it needs to generate a specialized version of the functions.

Let's see more in the details all the internal fields of the QuantumObject type:

julia
fieldnames(obj_type)
(:data, :type, :dimensions)
julia
σx_2.data
8×8 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries:
     ⋅          ⋅      1.0+0.0im  …      ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅    
 1.0+0.0im      ⋅          ⋅             ⋅          ⋅          ⋅    
     ⋅      1.0+0.0im      ⋅             ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅             ⋅      1.0+0.0im      ⋅    
     ⋅          ⋅          ⋅      …      ⋅          ⋅      1.0+0.0im
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅         1.0+0.0im      ⋅          ⋅
julia
σx_2.type
Operator

Operator is a synonym for OperatorQuantumObject.

julia
σx_2.dims
3-element SVector{3, Int64} with indices SOneTo(3):
 2
 2
 2

The dims field contains the dimensions of the subsystems (in this case, three subsystems with dimension 2 each). We can see that the type of dims is SVector instead of Vector. As we mentioned before, this is very useful in functions like ptrace. Let's do a simple example of reshaping an operator internally generated from some dims input:

julia
function reshape_operator_data(dims)
    op = Qobj(randn(prod(dims), prod(dims)), type=Operator, dims=dims)
    op_dims = op.dims
    op_data = op.data
    return reshape(op_data, vcat(op_dims, op_dims)...)
end

typeof(reshape_operator_data([2, 2, 2]))
Array{Float64, 6}

Which returns a tensor of size 2x2x2x2x2x2. Let's check the @code_warntype:

julia
@code_warntype reshape_operator_data([2, 2, 2])
MethodInstance for Main.var"Main".reshape_operator_data(::Vector{Int64})
  from reshape_operator_data(dims) @ Main.var"Main" type_stability.md:186
Arguments
  #self#::Core.Const(Main.var"Main".reshape_operator_data)
  dims::Vector{Int64}
Locals
  op_data::Matrix{Float64}
  op_dims::ANY
  op::QUANTUMOBJECT{OPERATORQUANTUMOBJECT, DIMENSIONS{_A, NTUPLE{N, SPACE}}, MATRIX{FLOAT64}} WHERE {_A, N}
Body::ANY
1 ─ %1  = Main.var"Main".randn::Core.Const(randn)
│   %2  = Main.var"Main".prod::Core.Const(prod)
│   %3  = (%2)(dims)::Int64
│   %4  = Main.var"Main".prod::Core.Const(prod)
│   %5  = (%4)(dims)::Int64
│   %6  = (%1)(%3, %5)::Matrix{Float64}
│   %7  = (:type, :dims)::Core.Const((:type, :dims))
│   %8  = Core.apply_type(Core.NamedTuple, %7)::Core.Const(NamedTuple{(:type, :dims)})
│   %9  = Main.var"Main".Operator::Core.Const(Operator)
│   %10 = Core.tuple(%9, dims)::Tuple{OperatorQuantumObject, Vector{Int64}}
│   %11 = (%8)(%10)::@NamedTuple{type::OperatorQuantumObject, dims::Vector{Int64}}
│   %12 = Main.var"Main".Qobj::Core.Const(QuantumObject)
│         (op = Core.kwcall(%11, %12, %6))
│   %14 = op::QUANTUMOBJECT{OPERATORQUANTUMOBJECT, DIMENSIONS{_A, NTUPLE{N, SPACE}}, MATRIX{FLOAT64}} WHERE {_A, N}
│         (op_dims = Base.getproperty(%14, :dims))
│   %16 = op::QUANTUMOBJECT{OPERATORQUANTUMOBJECT, DIMENSIONS{_A, NTUPLE{N, SPACE}}, MATRIX{FLOAT64}} WHERE {_A, N}
│         (op_data = Base.getproperty(%16, :data))
│   %18 = Main.var"Main".reshape::Core.Const(reshape)
│   %19 = op_data::Matrix{Float64}
│   %20 = Core.tuple(%19)::Tuple{Matrix{Float64}}
│   %21 = Main.var"Main".vcat::Core.Const(vcat)
│   %22 = op_dims::ANY
│   %23 = op_dims::ANY
│   %24 = (%21)(%22, %23)::ANY
│   %25 = Core._apply_iterate(Base.iterate, %18, %20, %24)::ANY
└──       return %25

We got a Any type, because the compiler doesn't know the size of the dims vector. We can fix this by using a Tuple (or SVector):

julia
typeof(reshape_operator_data((2, 2, 2)))
Array{Float64, 6}
julia
@code_warntype reshape_operator_data((2, 2, 2))
MethodInstance for Main.var"Main".reshape_operator_data(::Tuple{Int64, Int64, Int64})
  from reshape_operator_data(dims) @ Main.var"Main" type_stability.md:186
Arguments
  #self#::Core.Const(Main.var"Main".reshape_operator_data)
  dims::Tuple{Int64, Int64, Int64}
Locals
  op_data::Matrix{Float64}
  op_dims::SVector{3, Int64}
  op::QuantumObject{OperatorQuantumObject, Dimensions{3, Tuple{Space, Space, Space}}, Matrix{Float64}}
Body::Array{Float64, 6}
1 ─ %1  = Main.var"Main".randn::Core.Const(randn)
│   %2  = Main.var"Main".prod::Core.Const(prod)
│   %3  = (%2)(dims)::Int64
│   %4  = Main.var"Main".prod::Core.Const(prod)
│   %5  = (%4)(dims)::Int64
│   %6  = (%1)(%3, %5)::Matrix{Float64}
│   %7  = (:type, :dims)::Core.Const((:type, :dims))
│   %8  = Core.apply_type(Core.NamedTuple, %7)::Core.Const(NamedTuple{(:type, :dims)})
│   %9  = Main.var"Main".Operator::Core.Const(Operator)
│   %10 = Core.tuple(%9, dims)::Tuple{OperatorQuantumObject, Tuple{Int64, Int64, Int64}}
│   %11 = (%8)(%10)::@NamedTuple{type::OperatorQuantumObject, dims::Tuple{Int64, Int64, Int64}}
│   %12 = Main.var"Main".Qobj::Core.Const(QuantumObject)
│         (op = Core.kwcall(%11, %12, %6))
│   %14 = op::QuantumObject{OperatorQuantumObject, Dimensions{3, Tuple{Space, Space, Space}}, Matrix{Float64}}
│         (op_dims = Base.getproperty(%14, :dims))
│   %16 = op::QuantumObject{OperatorQuantumObject, Dimensions{3, Tuple{Space, Space, Space}}, Matrix{Float64}}
│         (op_data = Base.getproperty(%16, :data))
│   %18 = Main.var"Main".reshape::Core.Const(reshape)
│   %19 = op_data::Matrix{Float64}
│   %20 = Core.tuple(%19)::Tuple{Matrix{Float64}}
│   %21 = Main.var"Main".vcat::Core.Const(vcat)
│   %22 = op_dims::SVector{3, Int64}
│   %23 = op_dims::SVector{3, Int64}
│   %24 = (%21)(%22, %23)::SVector{6, Int64}
│   %25 = Core._apply_iterate(Base.iterate, %18, %20, %24)::Array{Float64, 6}
└──       return %25

Finally, let's look at the benchmarks

julia
@benchmark reshape_operator_data($[2, 2, 2])
BenchmarkTools.Trial: 10000 samples with 9 evaluations per sample.
 Range (min … max):  2.062 μs …  22.747 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.171 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.284 μs ± 912.686 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▄▇██▇▇▆▅▄▃▃▂▂▁      ▁▂▁▁                                   ▂
  ▇███████████████▇▆▇▆█████▇▇▅▆█████▇▄▄▆▆▆▆▆▆▅▂▄▄▂▅▅▇█▇▇▇▇▆▃▆ █
  2.06 μs      Histogram: log(frequency) by time      3.35 μs <

 Memory estimate: 1.50 KiB, allocs estimate: 29.
julia
@benchmark reshape_operator_data($((2, 2, 2)))
BenchmarkTools.Trial: 10000 samples with 580 evaluations per sample.
 Range (min … max):  206.867 ns … 86.539 μs  ┊ GC (min … max):  0.00% … 99.46%
 Time  (median):     221.100 ns              ┊ GC (median):     0.00%
 Time  (mean ± σ):   350.445 ns ±  1.325 μs  ┊ GC (mean ± σ):  15.66% ±  8.55%

   █▁                                                           
  ▃██▃▄▄▃▂▂▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▁▂▂▂▂▁▁▁▂▁▁▁▁▁▁▂▂▂▂▂▃▅▅▄▃▃▂▂▂▂▂▂▂▂▂▂ ▂
  207 ns          Histogram: frequency by time          540 ns <

 Memory estimate: 672 bytes, allocs estimate: 3.

Which is an innocuous but huge difference in terms of performance. Hence, we highly recommend using Tuple or SVector when defining the dimensions of a user-defined QuantumObject.

The use of Val in some QuantumToolbox.jl functions

In some functions of QuantumToolbox.jl, you may find the use of the Val type in the arguments. This is a trick to pass a value at compile time, and it is very useful to avoid type instabilities. Let's make a very simple example, where we want to create a Fock state |j of a given dimension N, and we give the possibility to create it as a sparse or dense vector. At first, we can write the function without using Val:

julia
function my_fock(N::Int, j::Int = 0; sparse::Bool = false)
    if sparse
        array = sparsevec([j + 1], [1.0 + 0im], N)
    else
        array = zeros(ComplexF64, N)
        array[j+1] = 1
    end
    return QuantumObject(array; type = Ket)
end
@show my_fock(2, 1)
@show my_fock(2, 1; sparse = true)
my_fock(2, 1) =
Quantum Object:   type=Ket   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
my_fock(2, 1; sparse = true) =
Quantum Object:   type=Ket   dims=[2]   size=(2,)
2-element SparseVector{ComplexF64, Int64} with 1 stored entry:
  [2]  =  1.0+0.0im

But it is immediately clear that the return type of this function is not clear, because it depends on the value of the sparse argument. We can check it using the @code_warntype macro:

julia
@code_warntype my_fock(2, 1)
MethodInstance for Main.var"Main".my_fock(::Int64, ::Int64)
  from my_fock(N::Int64, j::Int64; sparse) @ Main.var"Main" type_stability.md:229
Arguments
  #self#::Core.Const(Main.var"Main".my_fock)
  N::Int64
  j::Int64
Body::UNION{QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, VECTOR{COMPLEXF64}}, QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, SPARSEVECTOR{COMPLEXF64, INT64}}}
1 ─ %1 = Main.var"Main".:(var"#my_fock#1")::Core.Const(Main.var"Main".var"#my_fock#1")
│   %2 = (%1)(false, #self#, N, j)::UNION{QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, VECTOR{COMPLEXF64}}, QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, SPARSEVECTOR{COMPLEXF64, INT64}}}
└──      return %2
julia
@code_warntype my_fock(2, 1; sparse = true)
MethodInstance for Core.kwcall(::@NamedTuple{sparse::Bool}, ::typeof(Main.var"Main".my_fock), ::Int64, ::Int64)
  from kwcall(::NamedTuple, ::typeof(Main.var"Main".my_fock), N::Int64, j::Int64) @ Main.var"Main" type_stability.md:229
Arguments
  #s5::Core.Const(Core.kwcall)
  @_2::@NamedTuple{sparse::Bool}
  @_3::Core.Const(Main.var"Main".my_fock)
  N::Int64
  j::Int64
Locals
  sparse::Union{}
  @_7::Bool
Body::UNION{QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, VECTOR{COMPLEXF64}}, QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, SPARSEVECTOR{COMPLEXF64, INT64}}}
1 ──       Core.NewvarNode(:(sparse))
│          Core.NewvarNode(:(@_7))
│    %3  = Core.isdefined(@_2, :sparse)::Core.Const(true)
└───       goto #6 if not %3
2 ── %5  = Core.getfield(@_2, :sparse)::Bool
│    %6  = Main.var"Main".Bool::Core.Const(Bool)
│    %7  = (%5 isa %6)::Core.Const(true)
└───       goto #4 if not %7
3 ──       goto #5
4 ──       Core.Const(:(Main.var"Main".Bool))
│          Core.Const(:(%new(Core.TypeError, Symbol("keyword argument"), :sparse, %10, %5)))
└───       Core.Const(:(Core.throw(%11)))
5 ┄─       (@_7 = %5)
└───       goto #7
6 ──       Core.Const(:(@_7 = false))
7 ┄─ %16 = @_7::Bool
│    %17 = Base.keys(@_2)::Core.Const((:sparse,))
│    %18 = (:sparse,)::Core.Const((:sparse,))
│    %19 = Base.diff_names(%17, %18)::Core.Const(())
│    %20 = Base.isempty(%19)::Core.Const(true)
└───       goto #9 if not %20
8 ──       goto #10
9 ──       Core.Const(:(Base.kwerr(@_2, @_3, N, j)))
10 ┄ %24 = Main.var"Main".:(var"#my_fock#1")::Core.Const(Main.var"Main".var"#my_fock#1")
│    %25 = (%24)(%16, @_3, N, j)::UNION{QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, VECTOR{COMPLEXF64}}, QUANTUMOBJECT{KETQUANTUMOBJECT, DIMENSIONS{1, TUPLE{SPACE}}, SPARSEVECTOR{COMPLEXF64, INT64}}}
└───       return %25

We can fix this by using the Val type, where we enable the multiple dispatch of the function:

julia
getVal(::Val{N}) where N = N
function my_fock_good(N::Int, j::Int = 0; sparse::Val = Val(false))
    if getVal(sparse)
        array = zeros(ComplexF64, N)
        array[j+1] = 1
    else
        array = sparsevec([j + 1], [1.0 + 0im], N)
    end
    return QuantumObject(array; type = Ket)
end
@show my_fock_good(2, 1)
@show my_fock_good(2, 1; sparse = Val(true))
my_fock_good(2, 1) =
Quantum Object:   type=Ket   dims=[2]   size=(2,)
2-element SparseVector{ComplexF64, Int64} with 1 stored entry:
  [2]  =  1.0+0.0im
my_fock_good(2, 1; sparse = Val(true)) =
Quantum Object:   type=Ket   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im

And now the return type of the function is clear:

julia
@code_warntype my_fock_good(2, 1)
MethodInstance for Main.var"Main".my_fock_good(::Int64, ::Int64)
  from my_fock_good(N::Int64, j::Int64; sparse) @ Main.var"Main" type_stability.md:257
Arguments
  #self#::Core.Const(Main.var"Main".my_fock_good)
  N::Int64
  j::Int64
Body::QuantumObject{KetQuantumObject, Dimensions{1, Tuple{Space}}, SparseVector{ComplexF64, Int64}}
1 ─ %1 = Main.var"Main".:(var"#my_fock_good#2")::Core.Const(Main.var"Main".var"#my_fock_good#2")
│   %2 = Main.var"Main".Val(false)::Core.Const(Val{false}())
│   %3 = (%1)(%2, #self#, N, j)::QuantumObject{KetQuantumObject, Dimensions{1, Tuple{Space}}, SparseVector{ComplexF64, Int64}}
└──      return %3
julia
@code_warntype my_fock_good(2, 1; sparse = Val(true))
MethodInstance for Core.kwcall(::@NamedTuple{sparse::Val{true}}, ::typeof(Main.var"Main".my_fock_good), ::Int64, ::Int64)
  from kwcall(::NamedTuple, ::typeof(Main.var"Main".my_fock_good), N::Int64, j::Int64) @ Main.var"Main" type_stability.md:257
Arguments
  #s5::Core.Const(Core.kwcall)
  @_2::Core.Const((sparse = Val{true}(),))
  @_3::Core.Const(Main.var"Main".my_fock_good)
  N::Int64
  j::Int64
Locals
  sparse::Union{}
  @_7::Val{true}
Body::QuantumObject{KetQuantumObject, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}
1 ──       Core.NewvarNode(:(sparse))
│          Core.NewvarNode(:(@_7))
│    %3  = Core.isdefined(@_2, :sparse)::Core.Const(true)
└───       goto #6 if not %3
2 ── %5  = Core.getfield(@_2, :sparse)::Core.Const(Val{true}())
│    %6  = (%5 isa Main.var"Main".Val)::Core.Const(true)
└───       goto #4 if not %6
3 ──       goto #5
4 ──       Core.Const(:(%new(Core.TypeError, Symbol("keyword argument"), :sparse, Main.var"Main".Val, %5)))
└───       Core.Const(:(Core.throw(%9)))
5 ┄─       (@_7 = %5)
└───       goto #7
6 ──       Core.Const(:(@_7 = Main.var"Main".Val(false)))
7 ┄─ %14 = @_7::Core.Const(Val{true}())
│    %15 = Base.keys(@_2)::Core.Const((:sparse,))
│    %16 = (:sparse,)::Core.Const((:sparse,))
│    %17 = Base.diff_names(%15, %16)::Core.Const(())
│    %18 = Base.isempty(%17)::Core.Const(true)
└───       goto #9 if not %18
8 ──       goto #10
9 ──       Core.Const(:(Base.kwerr(@_2, @_3, N, j)))
10 ┄ %22 = Main.var"Main".:(var"#my_fock_good#2")::Core.Const(Main.var"Main".var"#my_fock_good#2")
│    %23 = (%22)(%14, @_3, N, j)::QuantumObject{KetQuantumObject, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}
└───       return %23

This is exactly how the current fock function is implemented in QuantumToolbox.jl. There are many other functions that support this feature, and we highly recommend using it when necessary.

Conclusions

In this page, we have seen the importance of type stability in Julia, and how to write efficient code in the context of QuantumToolbox.jl. We have seen that the internal structure of the QuantumObject type is already optimized for the compiler, and we have seen some practical examples of how to write efficient code. We have seen that the use of Vector should be avoided when the elements don't have the same type, and that the use of Tuple or SVector is highly recommended when the size of the array is known at compile time. Finally, we have seen the use of Val to pass values at compile time, to avoid type instabilities in some functions. ```