Allow only the struct to avoid type instabilities

Author: albertomercurio

docs/src/getting_started/type_stability.md

@@ -183,7 +183,7 @@ The `dims` field contains the dimensions of the subsystems (in this case, three
 
 ```@example type-stability
 function reshape_operator_data(dims)
-    op = Qobj(randn(prod(dims), prod(dims)), type=Operator, dims=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)...)
@@ -234,7 +234,7 @@ function my_fock(N::Int, j::Int = 0; sparse::Bool = false)
         array = zeros(ComplexF64, N)
         array[j+1] = 1
     end
-    return QuantumObject(array; type = Ket)
+    return QuantumObject(array; type = Ket())
 end
 @show my_fock(2, 1)
 @show my_fock(2, 1; sparse = true)
@@ -262,7 +262,7 @@ function my_fock_good(N::Int, j::Int = 0; sparse::Val = Val(false))
     else
         array = sparsevec([j + 1], [1.0 + 0im], N)
     end
-    return QuantumObject(array; type = Ket)
+    return QuantumObject(array; type = Ket())
 end
 @show my_fock_good(2, 1)
 @show my_fock_good(2, 1; sparse = Val(true))

docs/src/users_guide/QuantumObject/QuantumObject.md

@@ -54,7 +54,7 @@ Qobj(M, dims = SVector(2, 2)) # dims as StaticArrays.SVector (recommended)
     Please note that here we put the `dims` as a tuple `(2, 2)`. Although it supports also `Vector` type (`dims = [2, 2]`), it is recommended to use `Tuple` or `SVector` from [`StaticArrays.jl`](https://github.com/JuliaArrays/StaticArrays.jl) to improve performance. For a brief explanation on the impact of the type of `dims`, see the Section [The Importance of Type-Stability](@ref doc:Type-Stability).
 
 ```@example Qobj
-Qobj(rand(4, 4), type = SuperOperator)
+Qobj(rand(4, 4), type = SuperOperator())
 ```
 
 !!! note "Difference between `dims` and `size`"

docs/src/users_guide/states_and_operators.md

@@ -389,7 +389,7 @@ B = Qobj(rand(ComplexF64, N, N))
 mat2vec(A * ρ * B) ≈ spre(A) * spost(B) * mat2vec(ρ) ≈ sprepost(A, B) * mat2vec(ρ)
 ```
 
-In addition, dynamical generators on this extended space, often called Liouvillian superoperators, can be created using the [`liouvillian`](@ref) function. Each of these takes a Hamiltonian along with a list of collapse operators, and returns a [`type=SuperOperator`](@ref SuperOperator) object that can be exponentiated to find the superoperator for that evolution.
+In addition, dynamical generators on this extended space, often called Liouvillian superoperators, can be created using the [`liouvillian`](@ref) function. Each of these takes a Hamiltonian along with a list of collapse operators, and returns a [`QuantumObject`](@ref) of type [`SuperOperator`](@ref) that can be exponentiated to find the superoperator for that evolution.
 
 ```@example states_and_operators
 H = 10 * sigmaz()

src/negativity.jl

@@ -97,7 +97,7 @@ function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     ]
     return QuantumObject(
         reshape(permutedims(reshape(ρ.data, (dims..., dims...)), pt_idx), size(ρ)),
-        Operator,
+        Operator(),
         Dimensions(ρ.dimensions.to),
     )
 end
@@ -144,5 +144,5 @@ function _partial_transpose(
         end
     end
 
-    return QuantumObject(sparse(I_pt, J_pt, V_pt, M, N), Operator, ρ.dimensions)
+    return QuantumObject(sparse(I_pt, J_pt, V_pt, M, N), Operator(), ρ.dimensions)
 end

src/qobj/arithmetic_and_attributes.jl

@@ -69,7 +69,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                 )
                     check_dimensions(A, B)
                     QType = promote_op_type(A, B)
-                    return QType(A.data * B.data, Operator, A.dimensions)
+                    return QType(A.data * B.data, Operator(), A.dimensions)
                 end
             end
         else
@@ -82,7 +82,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                     QType = promote_op_type(A, B)
                     return QType(
                         A.data * B.data,
-                        Operator,
+                        Operator(),
                         GeneralDimensions(get_dimensions_to(A), get_dimensions_from(B)),
                     )
                 end
@@ -93,35 +93,35 @@ end
 
 function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:Dimensions})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
-    return QuantumObject(A.data * B.data, Ket, Dimensions(get_dimensions_to(A)))
+    return QuantumObject(A.data * B.data, Ket(), Dimensions(get_dimensions_to(A)))
 end
 function Base.:(*)(A::QuantumObject{Bra,<:Dimensions}, B::AbstractQuantumObject{Operator})
     check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
-    return QuantumObject(A.data * B.data, Bra, Dimensions(get_dimensions_from(B)))
+    return QuantumObject(A.data * B.data, Bra(), Dimensions(get_dimensions_from(B)))
 end
 function Base.:(*)(A::QuantumObject{Ket}, B::QuantumObject{Bra})
     check_dimensions(A, B)
-    return QuantumObject(A.data * B.data, Operator, A.dimensions) # to align with QuTiP, don't use kron(A, B) to do it.
+    return QuantumObject(A.data * B.data, Operator(), A.dimensions) # to align with QuTiP, don't use kron(A, B) to do it.
 end
 function Base.:(*)(A::QuantumObject{Bra}, B::QuantumObject{Ket})
     check_dimensions(A, B)
     return A.data * B.data
 end
 function Base.:(*)(A::AbstractQuantumObject{SuperOperator}, B::QuantumObject{Operator})
     check_dimensions(A, B)
-    return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dimensions)
+    return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator(), A.dimensions)
 end
 function Base.:(*)(A::QuantumObject{OperatorBra}, B::QuantumObject{OperatorKet})
     check_dimensions(A, B)
     return A.data * B.data
 end
 function Base.:(*)(A::AbstractQuantumObject{SuperOperator}, B::QuantumObject{OperatorKet})
     check_dimensions(A, B)
-    return QuantumObject(A.data * B.data, OperatorKet, A.dimensions)
+    return QuantumObject(A.data * B.data, OperatorKet(), A.dimensions)
 end
 function Base.:(*)(A::QuantumObject{OperatorBra}, B::AbstractQuantumObject{SuperOperator})
     check_dimensions(A, B)
-    return QuantumObject(A.data * B.data, OperatorBra, A.dimensions)
+    return QuantumObject(A.data * B.data, OperatorBra(), A.dimensions)
 end
 
 Base.:(^)(A::QuantumObject, n::T) where {T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
@@ -213,10 +213,10 @@ Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 """
 Base.adjoint(A::AbstractQuantumObject{OpType}) where {OpType<:Union{Operator,SuperOperator}} =
     get_typename_wrapper(A)(adjoint(A.data), A.type, adjoint(A.dimensions))
-Base.adjoint(A::QuantumObject{Ket}) = QuantumObject(adjoint(A.data), Bra, adjoint(A.dimensions))
-Base.adjoint(A::QuantumObject{Bra}) = QuantumObject(adjoint(A.data), Ket, adjoint(A.dimensions))
-Base.adjoint(A::QuantumObject{OperatorKet}) = QuantumObject(adjoint(A.data), OperatorBra, adjoint(A.dimensions))
-Base.adjoint(A::QuantumObject{OperatorBra}) = QuantumObject(adjoint(A.data), OperatorKet, adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{Ket}) = QuantumObject(adjoint(A.data), Bra(), adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{Bra}) = QuantumObject(adjoint(A.data), Ket(), adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{OperatorKet}) = QuantumObject(adjoint(A.data), OperatorBra(), adjoint(A.dimensions))
+Base.adjoint(A::QuantumObject{OperatorBra}) = QuantumObject(adjoint(A.data), OperatorKet(), adjoint(A.dimensions))
 
 @doc raw"""
     inv(A::AbstractQuantumObject)
@@ -521,7 +521,7 @@ function ptrace(QO::QuantumObject{Ket}, sel::Union{AbstractVector{Int},Tuple})
 
     _sort_sel = sort(SVector{length(sel),Int}(sel))
     ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator, dims = Dimensions(dkeep))
+    return QuantumObject(ρtr, type = Operator(), dims = Dimensions(dkeep))
 end
 
 ptrace(QO::QuantumObject{Bra}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
@@ -549,7 +549,7 @@ function ptrace(QO::QuantumObject{Operator}, sel::Union{AbstractVector{Int},Tupl
     dims = dimensions_to_dims(get_dimensions_to(QO))
     _sort_sel = sort(SVector{length(sel),Int}(sel))
     ρtr, dkeep = _ptrace_oper(QO.data, dims, _sort_sel)
-    return QuantumObject(ρtr, type = Operator, dims = Dimensions(dkeep))
+    return QuantumObject(ρtr, type = Operator(), dims = Dimensions(dkeep))
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 

src/qobj/eigsolve.jl

@@ -91,9 +91,9 @@ Base.iterate(res::EigsolveResult) = (res.values, Val(:vector_list))
 Base.iterate(res::EigsolveResult{T1,T2,Nothing}, ::Val{:vector_list}) where {T1,T2} =
     ([res.vectors[:, k] for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult{T1,T2,Operator}, ::Val{:vector_list}) where {T1,T2} =
-    ([QuantumObject(res.vectors[:, k], Ket, res.dimensions) for k in 1:length(res.values)], Val(:vectors))
+    ([QuantumObject(res.vectors[:, k], Ket(), res.dimensions) for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult{T1,T2,SuperOperator}, ::Val{:vector_list}) where {T1,T2} =
-    ([QuantumObject(res.vectors[:, k], OperatorKet, res.dimensions) for k in 1:length(res.values)], Val(:vectors))
+    ([QuantumObject(res.vectors[:, k], OperatorKet(), res.dimensions) for k in 1:length(res.values)], Val(:vectors))
 Base.iterate(res::EigsolveResult, ::Val{:vectors}) = (res.vectors, Val(:done))
 Base.iterate(res::EigsolveResult, ::Val{:done}) = nothing
 
@@ -408,7 +408,7 @@ function eigsolve_al(
     L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
     prob = mesolveProblem(
         L_evo,
-        QuantumObject(ρ0, type = Operator, dims = H.dimensions),
+        QuantumObject(ρ0, type = Operator(), dims = H.dimensions),
         [zero(T), T];
         params = params,
         progress_bar = Val(false),

src/qobj/functions.jl

@@ -208,7 +208,7 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                     _lazy_tensor_warning(A.data, B.data)
                     return QType(
                         kron(A.data, B.data),
-                        Operator,
+                        Operator(),
                         GeneralDimensions(
                             (get_dimensions_to(A)..., get_dimensions_to(B)...),
                             (get_dimensions_from(A)..., get_dimensions_from(B)...),
@@ -230,7 +230,7 @@ for AOpType in (:Ket, :Bra, :Operator)
                     _lazy_tensor_warning(A.data, B.data)
                     return QType(
                         kron(A.data, B.data),
-                        Operator,
+                        Operator(),
                         GeneralDimensions(
                             (get_dimensions_to(A)..., get_dimensions_to(B)...),
                             (get_dimensions_from(A)..., get_dimensions_from(B)...),
@@ -267,7 +267,7 @@ Convert a quantum object from vector ([`OperatorKet`](@ref)-type) to matrix ([`O
 !!! note
     `vector_to_operator` is a synonym of `vec2mat`.
 """
-vec2mat(A::QuantumObject{OperatorKet}) = QuantumObject(vec2mat(A.data), Operator, A.dimensions)
+vec2mat(A::QuantumObject{OperatorKet}) = QuantumObject(vec2mat(A.data), Operator(), A.dimensions)
 
 @doc raw"""
     mat2vec(A::QuantumObject)
@@ -278,7 +278,7 @@ Convert a quantum object from matrix ([`Operator`](@ref)-type) to vector ([`Oper
 !!! note
     `operator_to_vector` is a synonym of `mat2vec`.
 """
-mat2vec(A::QuantumObject{Operator}) = QuantumObject(mat2vec(A.data), OperatorKet, A.dimensions)
+mat2vec(A::QuantumObject{Operator}) = QuantumObject(mat2vec(A.data), OperatorKet(), A.dimensions)
 
 @doc raw"""
     mat2vec(A::AbstractMatrix)

src/qobj/operators.jl

@@ -48,7 +48,7 @@ function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, :
     # Because inv(Λ) ⋅ R has real and strictly positive elements, Q · Λ is therefore Haar distributed.
     Λ = diag(R) # take the diagonal elements of R
     Λ ./= abs.(Λ) # rescaling the elements
-    return QuantumObject(to_dense(Q * Diagonal(Λ)); type = Operator, dims = dimensions)
+    return QuantumObject(to_dense(Q * Diagonal(Λ)); type = Operator(), dims = dimensions)
 end
 function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, ::Val{:exp})
     N = prod(dimensions)
@@ -57,7 +57,7 @@ function rand_unitary(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}, :
     Z = randn(ComplexF64, N, N)
 
     # generate Hermitian matrix
-    H = QuantumObject((Z + Z') / 2; type = Operator, dims = dimensions)
+    H = QuantumObject((Z + Z') / 2; type = Operator(), dims = dimensions)
 
     return to_dense(exp(-1.0im * H))
 end
@@ -99,7 +99,7 @@ julia> fock(20, 3)' * a * fock(20, 4)
 2.0 + 0.0im
 ```
 """
-destroy(N::Int) = QuantumObject(spdiagm(1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator, N)
+destroy(N::Int) = QuantumObject(spdiagm(1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator(), N)
 
 @doc raw"""
     create(N::Int)
@@ -125,7 +125,7 @@ julia> fock(20, 4)' * a_d * fock(20, 3)
 2.0 + 0.0im
 ```
 """
-create(N::Int) = QuantumObject(spdiagm(-1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator, N)
+create(N::Int) = QuantumObject(spdiagm(-1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator(), N)
 
 @doc raw"""
     displace(N::Int, α::Number)
@@ -166,7 +166,7 @@ Bosonic number operator with Hilbert space cutoff `N`.
 
 This operator is defined as ``\hat{N}=\hat{a}^\dagger \hat{a}``, where ``\hat{a}`` is the bosonic annihilation operator.
 """
-num(N::Int) = QuantumObject(spdiagm(0 => Array{ComplexF64}(0:(N-1))), Operator, N)
+num(N::Int) = QuantumObject(spdiagm(0 => Array{ComplexF64}(0:(N-1))), Operator(), N)
 
 @doc raw"""
     position(N::Int)
@@ -222,7 +222,7 @@ function phase(N::Int, ϕ0::Real = 0)
     N_list = collect(0:(N-1))
     ϕ = ϕ0 .+ (2 * π / N) .* N_list
     states = [exp.((1.0im * ϕ[m]) .* N_list) ./ sqrt(N) for m in 1:N]
-    return QuantumObject(sum([ϕ[m] * states[m] * states[m]' for m in 1:N]); type = Operator, dims = N)
+    return QuantumObject(sum([ϕ[m] * states[m] * states[m]' for m in 1:N]); type = Operator(), dims = N)
 end
 
 @doc raw"""
@@ -276,36 +276,36 @@ function jmat(j::Real, ::Val{:x})
         throw(ArgumentError("The spin quantum number (j) must be a non-negative integer or half-integer."))
 
     σ = _jm(j)
-    return QuantumObject((σ' + σ) / 2, Operator, Int(J))
+    return QuantumObject((σ' + σ) / 2, Operator(), Int(J))
 end
 function jmat(j::Real, ::Val{:y})
     J = 2 * j + 1
     ((floor(J) != J) || (j < 0)) &&
         throw(ArgumentError("The spin quantum number (j) must be a non-negative integer or half-integer."))
 
     σ = _jm(j)
-    return QuantumObject((σ' - σ) / 2im, Operator, Int(J))
+    return QuantumObject((σ' - σ) / 2im, Operator(), Int(J))
 end
 function jmat(j::Real, ::Val{:z})
     J = 2 * j + 1
     ((floor(J) != J) || (j < 0)) &&
         throw(ArgumentError("The spin quantum number (j) must be a non-negative integer or half-integer."))
 
-    return QuantumObject(_jz(j), Operator, Int(J))
+    return QuantumObject(_jz(j), Operator(), Int(J))
 end
 function jmat(j::Real, ::Val{:+})
     J = 2 * j + 1
     ((floor(J) != J) || (j < 0)) &&
         throw(ArgumentError("The spin quantum number (j) must be a non-negative integer or half-integer."))
 
-    return QuantumObject(adjoint(_jm(j)), Operator, Int(J))
+    return QuantumObject(adjoint(_jm(j)), Operator(), Int(J))
 end
 function jmat(j::Real, ::Val{:-})
     J = 2 * j + 1
     ((floor(J) != J) || (j < 0)) &&
         throw(ArgumentError("The spin quantum number (j) must be a non-negative integer or half-integer."))
 
-    return QuantumObject(_jm(j), Operator, Int(J))
+    return QuantumObject(_jm(j), Operator(), Int(J))
 end
 jmat(j::Real, ::Val{T}) where {T} = throw(ArgumentError("Invalid spin operator: $(T)"))
 
@@ -428,8 +428,7 @@ Note that `type` can only be either [`Operator`](@ref) or [`SuperOperator`](@ref
 !!! note
     `qeye` is a synonym of `eye`.
 """
-function eye(N::Int; type = Operator, dims = nothing)
-    type = _get_type(type)
+function eye(N::Int; type = Operator(), dims = nothing)
     if dims isa Nothing
         dims = isa(type, Operator) ? N : isqrt(N)
     end
@@ -486,7 +485,7 @@ function _Jordan_Wigner(::Val{N}, j::Int, op::QuantumObject{Operator}) where {N}
     S = 2^(N - j)
     I_tensor = sparse((1.0 + 0.0im) * LinearAlgebra.I, S, S)
 
-    return QuantumObject(kron(Z_tensor, op.data, I_tensor); type = Operator, dims = ntuple(i -> 2, Val(N)))
+    return QuantumObject(kron(Z_tensor, op.data, I_tensor); type = Operator(), dims = ntuple(i -> 2, Val(N)))
 end
 
 @doc raw"""
@@ -495,7 +494,7 @@ end
 Generates the projection operator ``\hat{O} = |i \rangle\langle j|`` with Hilbert space dimension `N`.
 """
 projection(N::Int, i::Int, j::Int) =
-    QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type = Operator, dims = N)
+    QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type = Operator(), dims = N)
 
 @doc raw"""
     tunneling(N::Int, m::Int=1; sparse::Union{Bool,Val{<:Bool}}=Val(false))
@@ -518,9 +517,9 @@ function tunneling(N::Int, m::Int = 1; sparse::Union{Bool,Val} = Val(false))
 
     data = ones(ComplexF64, N - m)
     if getVal(sparse)
-        return QuantumObject(spdiagm(m => data, -m => data); type = Operator, dims = N)
+        return QuantumObject(spdiagm(m => data, -m => data); type = Operator(), dims = N)
     else
-        return QuantumObject(diagm(m => data, -m => data); type = Operator, dims = N)
+        return QuantumObject(diagm(m => data, -m => data); type = Operator(), dims = N)
     end
 end
 
@@ -551,9 +550,9 @@ where ``\omega = \exp(\frac{2 \pi i}{N})``.
 !!! warning "Beware of type-stability!"
     It is highly recommended to use `qft(dimensions)` with `dimensions` as `Tuple` or `SVector` from [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) to keep type stability. See the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-qft(dimensions::Int) = QuantumObject(_qft_op(dimensions), Operator, dimensions)
+qft(dimensions::Int) = QuantumObject(_qft_op(dimensions), Operator(), dimensions)
 qft(dimensions::Union{Dimensions,AbstractVector{Int},Tuple}) =
-    QuantumObject(_qft_op(prod(dimensions)), Operator, dimensions)
+    QuantumObject(_qft_op(prod(dimensions)), Operator(), dimensions)
 function _qft_op(N::Int)
     ω = exp(2.0im * π / N)
     arr = 0:(N-1)