Add LiouvilleSpace and restructure all structure of Dimensions and Spaces

CHANGELOG.md

@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Add `spre` and `spost` methods for `ComposedOperator` and cache propagator in every time evolution solver. ([#596])
+
 ## [v0.39.0]
 Release date: 2025-11-17
 
@@ -374,3 +376,4 @@ Release date: 2024-11-13
 [#588]: https://github.com/qutip/QuantumToolbox.jl/issues/588
 [#589]: https://github.com/qutip/QuantumToolbox.jl/issues/589
 [#591]: https://github.com/qutip/QuantumToolbox.jl/issues/591
+[#596]: https://github.com/qutip/QuantumToolbox.jl/issues/596

src/QuantumToolbox.jl

@@ -50,6 +50,7 @@ import SciMLOperators:
     ScalarOperator,
     ScaledOperator,
     AddedOperator,
+    ComposedOperator,
     IdentityOperator,
     update_coefficients!,
     concretize

src/qobj/arithmetic_and_attributes.jl

@@ -38,6 +38,8 @@ for op in (:(+), :(-), :(*))
     @eval begin
         function Base.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
             check_dimensions(A, B)
+            A.type != B.type &&
+                throw(DomainError((A.type, B.type), "The quantum objects should have the same type to do $op."))
             QType = promote_op_type(A, B)
             return QType($(op)(A.data, B.data), A.type, A.dimensions)
         end
@@ -50,54 +52,24 @@ for op in (:(+), :(-), :(*))
     end
 end
 
-function check_mul_dimensions(from::NTuple{NA,AbstractSpace}, to::NTuple{NB,AbstractSpace}) where {NA,NB}
-    (from != to) && throw(
-        DimensionMismatch(
-            "The quantum object with (right) dims = $(dimensions_to_dims(from)) can not multiply a quantum object with (left) dims = $(dimensions_to_dims(to)) on the right-hand side.",
-        ),
+for type in (:Operator, :SuperOperator)
+    @eval function Base.:(*)(
+        A::AbstractQuantumObject{$type,ProductDimensions},
+        B::AbstractQuantumObject{$type,ProductDimensions},
     )
-    return nothing
-end
-
-for ADimType in (:Dimensions, :GeneralDimensions)
-    for BDimType in (:Dimensions, :GeneralDimensions)
-        if ADimType == BDimType == :Dimensions
-            @eval begin
-                function Base.:(*)(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    check_dimensions(A, B)
-                    QType = promote_op_type(A, B)
-                    return QType(A.data * B.data, Operator(), A.dimensions)
-                end
-            end
-        else
-            @eval begin
-                function Base.:(*)(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
-                    QType = promote_op_type(A, B)
-                    return QType(
-                        A.data * B.data,
-                        Operator(),
-                        GeneralDimensions(get_dimensions_to(A), get_dimensions_from(B)),
-                    )
-                end
-            end
-        end
+        check_mul_dimensions(get_dimensions_from(A), get_dimensions_to(B))
+        QType = promote_op_type(A, B)
+        return QType(A.data * B.data, $type(), ProductDimensions(get_dimensions_to(A), get_dimensions_from(B)))
     end
 end
 
-function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket,<:Dimensions})
+function Base.:(*)(A::AbstractQuantumObject{Operator}, B::QuantumObject{Ket})
     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}, 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)
@@ -124,6 +96,15 @@ function Base.:(*)(A::QuantumObject{OperatorBra}, B::AbstractQuantumObject{Super
     return QuantumObject(A.data * B.data, OperatorBra(), A.dimensions)
 end
 
+function check_mul_dimensions(from::NTuple{NA,AbstractSpace}, to::NTuple{NB,AbstractSpace}) where {NA,NB}
+    (from != to) && throw(
+        DimensionMismatch(
+            "The quantum object with (right) dims = $(dimensions_to_dims(from)) can not multiply a quantum object with (left) dims = $(dimensions_to_dims(to)) on the right-hand side.",
+        ),
+    )
+    return nothing
+end
+
 Base.:(^)(A::QuantumObject, n::T) where {T<:Number} = QuantumObject(^(A.data, n), A.type, A.dimensions)
 Base.:(/)(A::AbstractQuantumObject, n::T) where {T<:Number} = get_typename_wrapper(A)(A.data / n, A.type, A.dimensions)
 
@@ -672,15 +653,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)
 
@@ -748,6 +729,5 @@ _dims_and_perm(::Operator, dims::SVector{N,Int}, order::AbstractVector{Int}, L::
 _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], isnothing(dimensions.from) ? nothing : dimensions.from[order])

src/qobj/dimensions.jl

@@ -2,103 +2,100 @@
 This file defines the Dimensions structures, which can describe composite Hilbert spaces.
 =#
 
-export AbstractDimensions, Dimensions, GeneralDimensions
+export AbstractDimensions, ProductDimensions
 
 abstract type AbstractDimensions{M,N} end
 
 @doc raw"""
-    struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N, N}
-        to::T
-    end
-
-A structure that describes the Hilbert [`Space`](@ref) of each subsystems.
-"""
-struct Dimensions{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)
-end
-function Dimensions(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)))
-end
-Dimensions(dims::Int) = Dimensions(Space(dims))
-Dimensions(dims::DimType) where {DimType<:AbstractSpace} = Dimensions((dims,))
-Dimensions(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 ProductDimensions{M,N,T1<:Tuple,T2<:Union{<:Tuple, Nothing}} <: AbstractDimensions{M,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 describes the left-hand side (`to`) and right-hand side (`from`) dimensions of a quantum object.
+
+The `from` field can be different from `nothing` only for non-square [`Operator`](@ref) and [`SuperOperator`](@ref) quantum objects.
 """
-struct GeneralDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
+struct ProductDimensions{M,N,T1<:Tuple,T2<:Union{<:Tuple,Nothing}} <: 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} =
-        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}
-    (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
+    function ProductDimensions(to::Union{AbstractVector,NTuple}, from::Union{NTuple,Nothing})
+        M = length(to)
+        N = isnothing(from) ? M : length(from)
 
-    _non_static_array_warning("dims[1]", dims[1])
-    _non_static_array_warning("dims[2]", dims[2])
+        _non_static_array_warning("dims", to)
+        isnothing(from) || _non_static_array_warning("dims", from)
 
-    L1 = length(dims[1])
-    L2 = length(dims[2])
-    (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."))
+        to_space = _dims_tuple_of_space(to)
+        from_space = _dims_tuple_of_space(from)
 
-    return GeneralDimensions(Tuple(Space.(dims[1])), Tuple(Space.(dims[2])))
+        new{M,N,typeof(to_space),typeof(from_space)}(to_space, from_space)
+    end
 end
+function ProductDimensions(dims::Union{AbstractVector,Tuple})
+    (length(dims) != 2) && throw(ArgumentError("Invalid dims = $dims"))
 
-_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{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N} =
-    GeneralDimensions(dims)
-_gen_dimensions(dims::Any) = Dimensions(dims)
+    return ProductDimensions(dims[1], dims[2])
+end
+
+ProductDimensions(dims::Union{Int,AbstractSpace}) = ProductDimensions((dims,), nothing)
+ProductDimensions(dims::Union{AbstractVector{<:Integer},NTuple{N,Integer}}) where {N} = ProductDimensions(dims, nothing)
+ProductDimensions(dims::Union{AbstractVector{<:AbstractSpace},NTuple{N,AbstractSpace}}) where {N} = ProductDimensions(dims, nothing)
+ProductDimensions(dims::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) =
-    SVector{2}(dimensions_to_dims(dimensions.to), dimensions_to_dims(dimensions.from))
-
+function dimensions_to_dims(dimensions::ProductDimensions)
+    dims_to = dimensions_to_dims(dimensions.to)
+    isnothing(dimensions.from) && return dims_to
+    dims_from = dimensions_to_dims(dimensions.from)
+    return SVector{2}(dims_to, dims_from)
+end
 dimensions_to_dims(::Nothing) = nothing # for EigsolveResult.dimensions = nothing
 
-Base.length(::AbstractDimensions{N}) where {N} = N
+hilbert_dimensions_to_size(dimensions::ProductDimensions) =
+    (hilbert_dimensions_to_size(dimensions.to), hilbert_dimensions_to_size(dimensions.from))
+hilbert_dimensions_to_size(dimensions::NTuple{N,AbstractSpace}) where {N} = prod(hilbert_dimensions_to_size, dimensions)
+hilbert_dimensions_to_size(dim::Int) = dim
+hilbert_dimensions_to_size(dimensions::Union{AbstractVector, NTuple{N,Integer}}) where {N} = prod(dimensions)
+hilbert_dimensions_to_size(::Nothing) = nothing
+
+liouvillian_dimensions_to_size(dimensions::ProductDimensions) =
+    (liouvillian_dimensions_to_size(dimensions.to), liouvillian_dimensions_to_size(dimensions.from))
+liouvillian_dimensions_to_size(dimensions::NTuple{N,AbstractSpace}) where {N} =
+    prod(liouvillian_dimensions_to_size, dimensions)
+liouvillian_dimensions_to_size(::Nothing) = nothing
 
-# 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)
+Base.length(::AbstractDimensions{N}) where {N} = N
 
-Base.transpose(dimensions::Dimensions) = dimensions
-Base.transpose(dimensions::GeneralDimensions) = GeneralDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions) = ProductDimensions(dimensions.from, dimensions.to) # switch `to` and `from`
+Base.transpose(dimensions::ProductDimensions{M,N,T1,Nothing}) where {M,N,T1<:Tuple} = dimensions
 Base.adjoint(dimensions::AbstractDimensions) = transpose(dimensions)
 
+Base.:(==)(dim1::ProductDimensions, dim2::ProductDimensions) = (dim1.to == dim2.to) && (dim1.from == dim2.from)
+
+_dims_tuple_of_space(dims::NTuple{N,AbstractSpace}) where {N} = dims
+function _dims_tuple_of_space(dims::Union{AbstractVector{<:Integer},SVector{M, <:Integer}, NTuple{M, Integer}}) where {M}
+    _non_static_array_warning("dims", dims)
+
+    N = length(dims)
+    N > 0 || throw(DomainError(N, "The length of `dims` must be larger or equal to 1."))
+
+    return ntuple(dim -> HilbertSpace(dims[dim]), Val(N))
+end
+_dims_tuple_of_space(::Nothing) = nothing
+
+_gen_dimensions(dims::AbstractDimensions) = dims
+_gen_dimensions(dims) = ProductDimensions(dims)
+
 # 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)
+function _get_dims_string(dimensions::ProductDimensions)
     dims = dimensions_to_dims(dimensions)
+    isnothing(dimensions.from) && return string(dims)
     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
+space_one_list(dimensions::NTuple{N,AbstractSpace}) where {N} =
+    ntuple(i -> one(dimensions[i]), Val(sum(length, dimensions)))

src/qobj/functions.jl

@@ -187,55 +187,34 @@ julia> a.dims, O.dims
 ```
 """
 function Base.kron(
-    A::AbstractQuantumObject{OpType,<:Dimensions},
-    B::AbstractQuantumObject{OpType,<:Dimensions},
-) where {OpType<:Union{Ket,Bra,Operator}}
+    A::AbstractQuantumObject{ObjType,<:ProductDimensions},
+    B::AbstractQuantumObject{ObjType,<:ProductDimensions},
+) where {ObjType<:Union{Ket, Bra, Operator}}
     QType = promote_op_type(A, B)
     _lazy_tensor_warning(A.data, B.data)
-    return QType(kron(A.data, B.data), A.type, Dimensions((A.dimensions.to..., B.dimensions.to...)))
-end
 
-# if A and B are both Operator but either one of them has GeneralDimensions
-for ADimType in (:Dimensions, :GeneralDimensions)
-    for BDimType in (:Dimensions, :GeneralDimensions)
-        if !(ADimType == BDimType == :Dimensions) # not for this case because it's already implemented
-            @eval begin
-                function Base.kron(
-                    A::AbstractQuantumObject{Operator,<:$ADimType},
-                    B::AbstractQuantumObject{Operator,<:$BDimType},
-                )
-                    QType = promote_op_type(A, B)
-                    _lazy_tensor_warning(A.data, B.data)
-                    return QType(
-                        kron(A.data, B.data),
-                        Operator(),
-                        GeneralDimensions(
-                            (get_dimensions_to(A)..., get_dimensions_to(B)...),
-                            (get_dimensions_from(A)..., get_dimensions_from(B)...),
-                        ),
-                    )
-                end
-            end
-        end
-    end
+    A_dims_from = get_dimensions_from(A)
+    B_dims_from = get_dimensions_from(B)
+    AB_dims_to = (get_dimensions_to(A)..., get_dimensions_to(B)...)
+    AB_dims_from = (A.dimensions.from === B.dimensions.from === nothing) ? nothing : (A_dims_from..., B_dims_from...)
+
+    return QType(kron(A.data, B.data), A.type, ProductDimensions(AB_dims_to, AB_dims_from))
 end
 
-# if A and B are different type (must return Operator with GeneralDimensions)
 for AOpType in (:Ket, :Bra, :Operator)
     for BOpType in (:Ket, :Bra, :Operator)
         if (AOpType != BOpType)
             @eval begin
                 function Base.kron(A::AbstractQuantumObject{$AOpType}, B::AbstractQuantumObject{$BOpType})
                     QType = promote_op_type(A, B)
                     _lazy_tensor_warning(A.data, B.data)
-                    return QType(
-                        kron(A.data, B.data),
-                        Operator(),
-                        GeneralDimensions(
-                            (get_dimensions_to(A)..., get_dimensions_to(B)...),
-                            (get_dimensions_from(A)..., get_dimensions_from(B)...),
-                        ),
-                    )
+
+                    A_dims_from = get_dimensions_from(A)
+                    B_dims_from = get_dimensions_from(B)
+                    AB_dims_to = (get_dimensions_to(A)..., get_dimensions_to(B)...)
+                    AB_dims_from = (A.dimensions.from === B.dimensions.from === nothing) ? nothing : (A_dims_from..., B_dims_from...)
+
+                    return QType(kron(A.data, B.data), Operator(), ProductDimensions(AB_dims_to, AB_dims_from))
                 end
             end
         end