Add LiouvilleSpace and restructure all structure of Dimensions and Spaces

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(hilbert_dimensions_to_size(dims)...), type=Operator(), dims=dims)
     op_dims = op.dims
     op_data = op.data
     return reshape(op_data, vcat(op_dims, op_dims)...)

docs/src/resources/api.md

@@ -17,10 +17,10 @@ end
 ## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
 
 ```@docs
-Space
+HilbertSpace
 EnrSpace
-Dimensions
-GeneralDimensions
+ProductDimensions
+GeneralProductDimensions
 AbstractQuantumObject
 Bra
 Ket
@@ -326,6 +326,8 @@ convert_unit
 row_major_reshape
 meshgrid
 enr_state_dictionaries
+hilbert_dimensions_to_size
+liouville_dimensions_to_size
 ```
 
 ## [Visualization](@id doc-API:Visualization)

ext/QuantumToolboxMakieExt.jl

@@ -226,7 +226,7 @@ function _plot_fock_distribution(
     kwargs...,
 ) where {SType<:Union{Bra,Ket,Operator}}
     ρ = ket2dm(ρ)
-    D = prod(ρ.dims)
+    D = hilbert_dimensions_to_size(ρ.dims)[1]
     isapprox(tr(ρ), 1, atol = 1e-4) || (@warn "The input ρ should be normalized.")
 
     xvec = 0:(D-1)

src/correlations.jl

@@ -113,7 +113,7 @@ function correlation_2op_2t(
     reverse::Bool = false,
     kwargs...,
 ) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator}}
-    C = eye(prod(H.dimensions), dims = H.dimensions)
+    C = eye(hilbert_dimensions_to_size(H.dimensions)[1], dims = H.dimensions)
     if reverse
         corr = correlation_3op_2t(H, ψ0, tlist, τlist, c_ops, A, B, C; kwargs...)
     else

src/entropy.jl

@@ -221,7 +221,7 @@ Calculate the [concurrence](https://en.wikipedia.org/wiki/Concurrence_(quantum_c
 - [Hill-Wootters1997](@citet)
 """
 function concurrence(ρ::QuantumObject{OpType}) where {OpType<:Union{Ket,Operator}}
-    (ρ.dimensions == Dimensions((Space(2), Space(2)))) || throw(
+    (ρ.dimensions == ProductDimensions((HilbertSpace(2), HilbertSpace(2)))) || throw(
         ArgumentError(
             "The `concurrence` only works for a two-qubit state, invalid dims = $(_get_dims_string(ρ.dimensions)).",
         ),

src/negativity.jl

@@ -76,7 +76,7 @@ function partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     (length(mask) != length(ρ.dimensions)) &&
         throw(ArgumentError("The length of \`mask\` should be equal to the length of \`ρ.dims\`."))
 
-    isa(ρ.dimensions, GeneralDimensions) &&
+    isa(ρ.dimensions, GeneralProductDimensions) &&
         (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
         throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
 
@@ -100,7 +100,7 @@ function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     return QuantumObject(
         reshape(permutedims(reshape(ρ.data, (dims_rev..., dims_rev...)), pt_idx), size(ρ)),
         Operator(),
-        Dimensions(ρ.dimensions.to),
+        ProductDimensions(ρ.dimensions.to),
     )
 end
 

src/qobj/arithmetic_and_attributes.jl

@@ -59,9 +59,9 @@ function check_mul_dimensions(from::NTuple{NA,AbstractSpace}, to::NTuple{NB,Abst
     return nothing
 end
 
-for ADimType in (:Dimensions, :GeneralDimensions)
-    for BDimType in (:Dimensions, :GeneralDimensions)
-        if ADimType == BDimType == :Dimensions
+for ADimType in (:ProductDimensions, :GeneralProductDimensions)
+    for BDimType in (:ProductDimensions, :GeneralProductDimensions)
+        if ADimType == BDimType == :ProductDimensions
             @eval begin
                 function Base.:(*)(
                     A::AbstractQuantumObject{Operator,<:$ADimType},
@@ -83,21 +83,21 @@ for ADimType in (:Dimensions, :GeneralDimensions)
                     return QType(
                         A.data * B.data,
                         Operator(),
-                        GeneralDimensions(get_dimensions_to(A), get_dimensions_from(B)),
+                        GeneralProductDimensions(get_dimensions_to(A), get_dimensions_from(B)),
                     )
                 end
             end
         end
     end
 end
 
-function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:Dimensions})
+function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:ProductDimensions})
     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(), ProductDimensions(get_dimensions_to(A)))
 end
-function Base.:(*)(A::QuantumObject{Bra,<:Dimensions}, B::AbstractQuantumObject{Operator})
+function Base.:(*)(A::QuantumObject{Bra,<:ProductDimensions}, 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(), ProductDimensions(get_dimensions_from(B)))
 end
 function Base.:(*)(A::QuantumObject{Ket}, B::QuantumObject{Bra})
     check_dimensions(A, B)
@@ -523,7 +523,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 = ProductDimensions(dkeep))
 end
 
 ptrace(QO::QuantumObject{Bra}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
@@ -532,7 +532,7 @@ function ptrace(QO::QuantumObject{Operator}, sel::Union{AbstractVector{Int},Tupl
     any(s -> s isa EnrSpace, QO.dimensions.to) && throw(ArgumentError("ptrace does not support EnrSpace"))
 
     # TODO: support for special cases when some of the subsystems have same `to` and `from` space
-    isa(QO.dimensions, GeneralDimensions) &&
+    isa(QO.dimensions, GeneralProductDimensions) &&
         (get_dimensions_to(QO) != get_dimensions_from(QO)) &&
         throw(ArgumentError("Invalid partial trace for dims = $(_get_dims_string(QO.dimensions))"))
 
@@ -553,7 +553,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 = ProductDimensions(dkeep))
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 
@@ -672,15 +672,15 @@ Get the coherence value ``\alpha`` by measuring the expectation value of the des
 It returns both ``\alpha`` and the corresponding state with the coherence removed: ``\ket{\delta_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \ket{\psi}`` for a pure state, and ``\hat{\rho_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \hat{\rho} \exp ( -\bar{\alpha} \hat{a} + \alpha \hat{a}^\dagger )`` for a density matrix. These states correspond to the quantum fluctuations around the coherent state ``\ket{\alpha}`` or ``|\alpha\rangle\langle\alpha|``.
 """
 function get_coherence(ψ::QuantumObject{Ket})
-    a = destroy(prod(ψ.dimensions))
+    a = destroy(hilbert_dimensions_to_size(ψ.dimensions)[1])
     α = expect(a, ψ)
     D = exp(α * a' - conj(α) * a)
 
     return α, D' * ψ
 end
 
 function get_coherence(ρ::QuantumObject{Operator})
-    a = destroy(prod(ρ.dimensions))
+    a = destroy(hilbert_dimensions_to_size(ρ.dimensions)[1])
     α = expect(a, ρ)
     D = exp(α * a' - conj(α) * a)
 
@@ -740,14 +740,14 @@ end
 _dims_and_perm(::ObjType, dims::SVector{N,Int}, order::AbstractVector{Int}, L::Int) where {ObjType<:Union{Ket,Bra},N} =
     reverse(dims), reverse((L + 1) .- order)
 
-# if dims originates from Dimensions
+# if dims originates from ProductDimensions
 _dims_and_perm(::Operator, dims::SVector{N,Int}, order::AbstractVector{Int}, L::Int) where {N} =
     reverse(vcat(dims, dims)), reverse((2 * L + 1) .- vcat(order, order .+ L))
 
-# if dims originates from GeneralDimensions
+# if dims originates from GeneralProductDimensions
 _dims_and_perm(::Operator, dims::SVector{2,SVector{N,Int}}, order::AbstractVector{Int}, L::Int) where {N} =
     reverse(vcat(dims[2], dims[1])), reverse((2 * L + 1) .- vcat(order, order .+ L))
 
-_order_dimensions(dimensions::Dimensions, order::AbstractVector{Int}) = Dimensions(dimensions.to[order])
-_order_dimensions(dimensions::GeneralDimensions, order::AbstractVector{Int}) =
-    GeneralDimensions(dimensions.to[order], dimensions.from[order])
+_order_dimensions(dimensions::ProductDimensions, order::AbstractVector{Int}) = ProductDimensions(dimensions.to[order])
+_order_dimensions(dimensions::GeneralProductDimensions, order::AbstractVector{Int}) =
+    GeneralProductDimensions(dimensions.to[order], dimensions.from[order])

src/qobj/dimensions.jl

@@ -1,56 +1,57 @@
 #=
-This file defines the Dimensions structures, which can describe composite Hilbert spaces.
+This file defines the ProductDimensions structures, which can describe composite Hilbert spaces.
 =#
 
-export AbstractDimensions, Dimensions, GeneralDimensions
+export AbstractDimensions, ProductDimensions, GeneralProductDimensions
+export hilbert_dimensions_to_size, liouville_dimensions_to_size
 
 abstract type AbstractDimensions{M,N} end
 
 @doc raw"""
-    struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
+    struct ProductDimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
         to::T
     end
 
-A structure that describes the Hilbert [`Space`](@ref) of each subsystems.
+A structure that embodies the [`AbstractSpace`](@ref) of each subsystem in a composite Hilbert space.
 """
-struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N,N}
+struct ProductDimensions{N,T<:Tuple} <: AbstractDimensions{N,N}
     to::T
 
     # make sure the elements in the tuple are all AbstractSpace
-    Dimensions(to::NTuple{N,AbstractSpace}) where {N} = new{N,typeof(to)}(to)
+    ProductDimensions(to::NTuple{N,AbstractSpace}) where {N} = new{N,typeof(to)}(to)
 end
-function Dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
+function ProductDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
     _non_static_array_warning("dims", dims)
     L = length(dims)
     (L > 0) || throw(DomainError(dims, "The argument dims must be of non-zero length"))
 
-    return Dimensions(Tuple(Space.(dims)))
+    return ProductDimensions(Tuple(HilbertSpace.(dims)))
 end
-Dimensions(dims::Int) = Dimensions(Space(dims))
-Dimensions(dims::DimType) where {DimType<:AbstractSpace} = Dimensions((dims,))
-Dimensions(dims::Any) = throw(
+ProductDimensions(dims::Int) = ProductDimensions(HilbertSpace(dims))
+ProductDimensions(dims::DimType) where {DimType<:AbstractSpace} = ProductDimensions((dims,))
+ProductDimensions(dims::Any) = throw(
     ArgumentError(
         "The argument dims must be a Tuple or a StaticVector of non-zero length and contain only positive integers.",
     ),
 )
 
 @doc raw"""
-    struct GeneralDimensions{N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{N}
+    struct GeneralProductDimensions{N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{N}
         to::T1
         from::T2
     end
 
-A structure that describes the left-hand side (`to`) and right-hand side (`from`) Hilbert [`Space`](@ref) of an [`Operator`](@ref).
+A structure that embodies the left-hand side (`to`) and right-hand side (`from`) [`AbstractSpace`](@ref) of a quantum object.
 """
-struct GeneralDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
+struct GeneralProductDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
     to::T1   # space acting on the left
     from::T2 # space acting on the right
 
     # make sure the elements in the tuple are all AbstractSpace
-    GeneralDimensions(to::NTuple{M,AbstractSpace}, from::NTuple{N,AbstractSpace}) where {M,N} =
+    GeneralProductDimensions(to::NTuple{M,AbstractSpace}, from::NTuple{N,AbstractSpace}) where {M,N} =
         new{M,N,typeof(to),typeof(from)}(to, from)
 end
-function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N}
+function GeneralProductDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N}
     (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
 
     _non_static_array_warning("dims[1]", dims[1])
@@ -61,44 +62,86 @@ function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<
     (L1 > 0) || throw(DomainError(L1, "The length of `dims[1]` must be larger or equal to 1."))
     (L2 > 0) || throw(DomainError(L2, "The length of `dims[2]` must be larger or equal to 1."))
 
-    return GeneralDimensions(Tuple(Space.(dims[1])), Tuple(Space.(dims[2])))
+    return GeneralProductDimensions(Tuple(HilbertSpace.(dims[1])), Tuple(HilbertSpace.(dims[2])))
 end
 
 _gen_dimensions(dims::AbstractDimensions) = dims
-_gen_dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N} = Dimensions(dims)
+_gen_dimensions(dims::Union{AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} = ProductDimensions(dims)
 _gen_dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N} =
-    GeneralDimensions(dims)
-_gen_dimensions(dims::Any) = Dimensions(dims)
+    GeneralProductDimensions(dims)
+_gen_dimensions(dims::Any) = ProductDimensions(dims)
 
 # obtain dims in the type of SVector with integers
 dimensions_to_dims(dimensions::NTuple{N,AbstractSpace}) where {N} = vcat(map(dimensions_to_dims, dimensions)...)
-dimensions_to_dims(dimensions::Dimensions) = dimensions_to_dims(dimensions.to)
-dimensions_to_dims(dimensions::GeneralDimensions) =
+dimensions_to_dims(dimensions::ProductDimensions) = dimensions_to_dims(dimensions.to)
+dimensions_to_dims(dimensions::GeneralProductDimensions) =
     SVector{2}(dimensions_to_dims(dimensions.to), dimensions_to_dims(dimensions.from))
 
 dimensions_to_dims(::Nothing) = nothing # for EigsolveResult.dimensions = nothing
 
 Base.length(::AbstractDimensions{N}) where {N} = N
 
-# need to specify return type `Int` for `_get_space_size`
-# otherwise the type of `prod(::Dimensions)` will be unstable
-_get_space_size(s::AbstractSpace)::Int = s.size
-Base.prod(dims::Dimensions) = prod(dims.to)
-Base.prod(spaces::NTuple{N,AbstractSpace}) where {N} = prod(_get_space_size, spaces)
+"""
+    hilbert_dimensions_to_size(dimensions)
+
+Returns the matrix dimensions `(m, n)` of an [`Operator`](@ref) with the given `dimensions`.
+
+For [`ProductDimensions`](@ref), returns `(m, m)` where `m` is the product of all subsystem Hilbert space dimensions.
+For [`GeneralProductDimensions`](@ref), returns `(m, n)` where `m` is the product of the `to` dimensions
+and `n` is the product of the `from` dimensions.
+
+If `dimensions` is an `Integer` or a vector/tuple of `Integer`s, it is automatically treated as [`ProductDimensions`](@ref).
+"""
+function hilbert_dimensions_to_size(dimensions::ProductDimensions)
+    m = prod(hilbert_dimensions_to_size, dimensions.to)
+    return (m, m)
+end
+function hilbert_dimensions_to_size(dimensions::GeneralProductDimensions)
+    m = prod(hilbert_dimensions_to_size, dimensions.to)
+    n = prod(hilbert_dimensions_to_size, dimensions.from)
+    return (m, n)
+end
+hilbert_dimensions_to_size(dimensions::Union{<:Integer,AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} =
+    hilbert_dimensions_to_size(ProductDimensions(dimensions))
+
+"""
+    liouville_dimensions_to_size(dimensions)
+
+Returns the matrix dimensions `(m, n)` of a [`SuperOperator`](@ref) with the given `dimensions`.
+
+For [`ProductDimensions`](@ref), returns `(m, m)` where `m` is the product of all subsystem Liouville space dimensions
+(each Hilbert dimension `d` contributes `d²` to the product).
+For [`GeneralProductDimensions`](@ref), returns `(m, n)` where `m` is the product of the `to` dimensions
+and `n` is the product of the `from` dimensions.
+
+If `dimensions` is an `Integer` or a vector/tuple of `Integer`s, it is automatically treated as [`ProductDimensions`](@ref).
+"""
+function liouville_dimensions_to_size(dimensions::ProductDimensions)
+    m = prod(liouville_dimensions_to_size, dimensions.to)
+    return (m, m)
+end
+function liouville_dimensions_to_size(dimensions::GeneralProductDimensions)
+    m = prod(liouville_dimensions_to_size, dimensions.to)
+    n = prod(liouville_dimensions_to_size, dimensions.from)
+    return (m, n)
+end
+liouville_dimensions_to_size(dimensions::Union{<:Integer,AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} =
+    liouville_dimensions_to_size(ProductDimensions(dimensions))
 
-Base.transpose(dimensions::Dimensions) = dimensions
-Base.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions) = dimensions
+Base.transpose(dimensions::GeneralProductDimensions) = GeneralProductDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
 Base.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
 
 # this is used to show `dims` for Qobj and QobjEvo
-_get_dims_string(dimensions::Dimensions) = string(dimensions_to_dims(dimensions))
-function _get_dims_string(dimensions::GeneralDimensions)
+_get_dims_string(dimensions::ProductDimensions) = string(dimensions_to_dims(dimensions))
+function _get_dims_string(dimensions::GeneralProductDimensions)
     dims = dimensions_to_dims(dimensions)
     return "[$(string(dims[1])), $(string(dims[2]))]"
 end
 _get_dims_string(::Nothing) = "nothing" # for EigsolveResult.dimensions = nothing
 
-Base.:(==)(dim1::Dimensions, dim2::Dimensions) = dim1.to == dim2.to
-Base.:(==)(dim1::GeneralDimensions, dim2::GeneralDimensions) = (dim1.to == dim2.to) && (dim1.from == dim2.from)
-Base.:(==)(dim1::Dimensions, dim2::GeneralDimensions) = false
-Base.:(==)(dim1::GeneralDimensions, dim2::Dimensions) = false
+Base.:(==)(dim1::ProductDimensions, dim2::ProductDimensions) = dim1.to == dim2.to
+Base.:(==)(dim1::GeneralProductDimensions, dim2::GeneralProductDimensions) =
+    (dim1.to == dim2.to) && (dim1.from == dim2.from)
+Base.:(==)(dim1::ProductDimensions, dim2::GeneralProductDimensions) = false
+Base.:(==)(dim1::GeneralProductDimensions, dim2::ProductDimensions) = false

src/qobj/eigsolve.jl

@@ -45,7 +45,7 @@ julia> λ
   1.0 + 0.0im
 
 julia> ψ
-2-element Vector{QuantumObject{Ket, Dimensions{1, Tuple{Space}}, Vector{ComplexF64}}}:
+2-element Vector{QuantumObject{Ket, ProductDimensions{1, Tuple{HilbertSpace}}, Vector{ComplexF64}}}:
 
 Quantum Object:   type=Ket()   dims=[2]   size=(2,)
 2-element Vector{ComplexF64}:
@@ -424,7 +424,7 @@ end
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         alg::AbstractODEAlgorithm = DP5(),
         params::NamedTuple = NamedTuple(),
-        ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
+        ρ0::AbstractMatrix = rand_dm(hilbert_dimensions_to_size(H.dimensions)[1]).data,
         eigvals::Int = 1,
         krylovdim::Int = min(10, size(H, 1)),
         maxiter::Int = 200,
@@ -467,7 +467,7 @@ function eigsolve_al(
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::AbstractODEAlgorithm = DP5(),
     params::NamedTuple = NamedTuple(),
-    ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
+    ρ0::AbstractMatrix = rand_dm(hilbert_dimensions_to_size(H.dimensions)[1]).data,
     eigvals::Int = 1,
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,