Type stability

[nmheim]

There are a few places with sub-optimal type definitions, like e.g.:

struct Shadow
    n_qubits::Int
    gate::Array{Yao.AbstractBlock}
    measure::Array{Int8}  # is this always a vector or can this be any dimensionality?
end

should really be

struct Shadow{G<:AbstractBlock}
  n_qubits::Int
  gate::G
  measure::Vector{Int8}
end

in order to make Shadow a proper concrete type that the julia compiler can work with optimally. without the parametric type G julia will always have to do dynamic dispatch (during runtime = slow) because the concrete type of the Shadow is not inferable (the gate can be any type of AbstractBlock)

[nmheim]

There is a similar problem with pre-allocating sets/arrays like below

QuantumMeasurements.jl/src/state_reconstruction.jl

Lines 54 to 56 in f985f5d

	ρ = Set{Matrix}()
	for z in keys(f)
	ρj = Array{Matrix}(undef, n)

which results in type unstable code, because Matrix is not a concrete type (as opposed to e.g. Matrix{Float64}. We should check if this is impacting performance, because julia has to runtime dispatch on functions things like this.
Things like this can be fixed by using functional primitives like map which will automatically infer the correct type for the resulting array.

[nmheim]

A more concrete example with a part of the function above. One can rewrite a for loop with manual preallocation to a call to map which will automatically infer the correct type for the resulting array:

ρj = Array{Matrix}(undef, n)
for j = 1:n
    vec = b[j] == 0 ? [1, 0] : [0, 1]

    #Compute the inverse of the channel measurement to return an estimation of the state's density matrix
    Π(A::Matrix) = tr(A) * [1 0; 0 1] / 2
    U = Matrix(mat(gate[j]))

    if z[j] == 0
        P = Π(U' * vec * vec' * U)
    else
        P = U' * vec * vec' * U - Π(U' * vec * vec' * U)
    end
    ρj[j] = P
end

# is equivalent to:
ρj = map(1:n) do j
    vec = b[j] == 0 ? [1, 0] : [0, 1]

    #Compute the inverse of the channel measurement to return an estimation of
    #the state's density matrix
    Π(A::Matrix) = tr(A) * [1 0; 0 1] / 2
    U = Matrix(mat(gate[j]))

    if z[j] == 0
        P = Π(U' * vec * vec' * U)
    else
        P = U' * vec * vec' * U - Π(U' * vec * vec' * U)
    end
end