Add tests for QObjEvo and time evolution

Author: albertomercurio

src/qobj/arithmetic_and_attributes.jl

@@ -38,7 +38,7 @@ for op in (:(+), :(-), :(*))
     @eval begin
         function LinearAlgebra.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
             check_dims(A, B)
-            QType = promote_type(A, B)
+            QType = promote_op_type(A, B)
             return QType($(op)(A.data, B.data), A.type, A.dims)
         end
         LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dims)

src/qobj/boolean_functions.jl

@@ -10,7 +10,7 @@ export isunitary, isconstant
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`BraQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isbra(A::QuantumObject) = A.type isa BraQuantumObject
+isbra(A::QuantumObject) = isbra(typeof(A))
 isbra(::Type{QuantumObject{DT,BraQuantumObject,N}}) where {DT,N} = true
 isbra(A) = false # default case
 
@@ -19,7 +19,7 @@ isbra(A) = false # default case
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`KetQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isket(A::QuantumObject) = A.type isa KetQuantumObject
+isket(A::QuantumObject) = isket(typeof(A))
 isket(::Type{QuantumObject{DT,KetQuantumObject,N}}) where {DT,N} = true
 isket(A) = false # default case
 
@@ -28,16 +28,16 @@ isket(A) = false # default case
 
 Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`OperatorQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoper(A::AbstractQuantumObject) = A.type isa OperatorQuantumObject
-isoper(::Type{QuantumObject{DT,OperatorQuantumObject,N}}) where {DT,N} = true
+isoper(A::AbstractQuantumObject) = isoper(typeof(A))
+isoper(::Type{<:AbstractQuantumObject{DT,OperatorQuantumObject,N}}) where {DT,N} = true
 isoper(A) = false # default case
 
 @doc raw"""
     isoperbra(A::QuantumObject)
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorBraQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoperbra(A::QuantumObject) = A.type isa OperatorBraQuantumObject
+isoperbra(A::QuantumObject) = isoperbra(typeof(A))
 isoperbra(::Type{QuantumObject{DT,OperatorBraQuantumObject,N}}) where {DT,N} = true
 isoperbra(A) = false # default case
 
@@ -46,7 +46,7 @@ isoperbra(A) = false # default case
 
 Checks if the [`QuantumObject`](@ref) `A` is a [`OperatorKetQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-isoperket(A::QuantumObject) = A.type isa OperatorKetQuantumObject
+isoperket(A::QuantumObject) = isoperket(typeof(A))
 isoperket(::Type{QuantumObject{DT,OperatorKetQuantumObject,N}}) where {DT,N} = true
 isoperket(A) = false # default case
 
@@ -55,8 +55,8 @@ isoperket(A) = false # default case
 
 Checks if the [`AbstractQuantumObject`](@ref) `A` is a [`SuperOperatorQuantumObject`](@ref). Default case returns `false` for any other inputs.
 """
-issuper(A::AbstractQuantumObject) = A.type isa SuperOperatorQuantumObject
-issuper(::Type{QuantumObject{DT,SuperOperatorQuantumObject,N}}) where {DT,N} = true
+issuper(A::AbstractQuantumObject) = issuper(typeof(A))
+issuper(::Type{<:AbstractQuantumObject{DT,SuperOperatorQuantumObject,N}}) where {DT,N} = true
 issuper(A) = false # default case
 
 @doc raw"""

src/qobj/functions.jl

@@ -193,7 +193,7 @@ function LinearAlgebra.kron(
     A::AbstractQuantumObject{DT1,OpType},
     B::AbstractQuantumObject{DT2,OpType},
 ) where {DT1,DT2,OpType<:Union{KetQuantumObject,BraQuantumObject,OperatorQuantumObject}}
-    QType = promote_type(A, B)
+    QType = promote_op_type(A, B)
     return QType(kron(A.data, B.data), A.type, vcat(A.dims, B.dims))
 end
 LinearAlgebra.kron(A::AbstractQuantumObject) = A

src/qobj/quantum_object_evo.jl

@@ -42,7 +42,7 @@ coef1 (generic function with 1 method)
 julia> coef2(p, t) = sin(t)
 coef2 (generic function with 1 method)
 
-julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
+julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
 Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
 (ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
 ```
@@ -67,7 +67,7 @@ coef1 (generic function with 1 method)
 julia> coef2(p, t) = sin(p.ω2 * t)
 coef2 (generic function with 1 method)
 
-julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
+julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
 Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
 (ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
 
@@ -270,15 +270,149 @@ function _make_SciMLOperator(op::QuantumObject, α; f::Function = identity)
     return MatrixOperator(α * f(op.data))
 end
 
-function (QO::QuantumObjectEvolution)(p, t)
+@doc raw"""
+    (A::QuantumObjectEvolution)(ψout, ψin, p, t)
+
+Apply the time-dependent [`QuantumObjectEvolution`](@ref) object `A` to the input state `ψin` at time `t` with parameters `p`. The output state is stored in `ψout`. This function mimics the behavior of a `AbstractSciMLOperator` object.
+
+# Arguments
+- `ψout::QuantumObject`: The output state. It must have the same type as `ψin`.
+- `ψin::QuantumObject`: The input state. It must be either a [`KetQuantumObject`](@ref) or a [`OperatorKetQuantumObject`](@ref).
+- `p`: The parameters of the time-dependent coefficients.
+- `t`: The time at which the coefficients are evaluated.
+
+# Returns
+- `ψout::QuantumObject`: The output state.
+
+# Examples
+```
+julia> a = destroy(20)
+Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
+⎡⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⎦
+
+julia> coef1(p, t) = sin(t)
+coef1 (generic function with 1 method)
+
+julia> coef2(p, t) = cos(t)
+coef2 (generic function with 1 method)
+
+julia> A = QobjEvo(((a, coef1), (a', coef2)))
+Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=true
+(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))
+
+julia> ψ1 = fock(20, 3)
+Quantum Object:   type=Ket   dims=[20]   size=(20,)
+20-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 1.0 + 0.0im
+ 0.0 + 0.0im
+     ⋮
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> ψ2 = zero_ket(20)
+Quantum Object:   type=Ket   dims=[20]   size=(20,)
+20-element Vector{ComplexF64}:
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+     ⋮
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+ 0.0 + 0.0im
+
+julia> A(ψ2, ψ1, nothing, 0.1)
+20-element Vector{ComplexF64}:
+                0.0 + 0.0im
+                0.0 + 0.0im
+ 0.1729165499254989 + 0.0im
+                0.0 + 0.0im
+ 1.9900083305560516 + 0.0im
+                    ⋮
+                0.0 + 0.0im
+                0.0 + 0.0im
+                0.0 + 0.0im
+                0.0 + 0.0im
+```
+"""
+function (A::QuantumObjectEvolution)(
+    ψout::QuantumObject{DT1,QobjType},
+    ψin::QuantumObject{DT2,QobjType},
+    p,
+    t,
+) where {DT1,DT2,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    check_dims(ψout, ψin)
+    check_dims(ψout, A)
+
+    if isoper(A) && isoperket(ψin)
+        throw(
+            ArgumentError(
+                "The input state must be a KetQuantumObject if the QuantumObjectEvolution object is an Operator.",
+            ),
+        )
+    elseif issuper(A) && isket(ψin)
+        throw(
+            ArgumentError(
+                "The input state must be an OperatorKetQuantumObject if the QuantumObjectEvolution object is a SuperOperator.",
+            ),
+        )
+    end
+
+    A.data(ψout.data, ψin.data, p, t)
+
+    return ψout
+end
+
+@doc raw"""
+    (A::QuantumObjectEvolution)(ψ, p, t)
+
+Apply the time-dependent [`QuantumObjectEvolution`](@ref) object `A` to the input state `ψ` at time `t` with parameters `p`. Out-of-place version of [`(A::QuantumObjectEvolution)(ψout, ψin, p, t)`](@ref). The output state is stored in a new [`QuantumObject`](@ref) object. This function mimics the behavior of a `AbstractSciMLOperator` object.
+"""
+function (A::QuantumObjectEvolution)(
+    ψ::QuantumObject{DT,QobjType},
+    p,
+    t,
+) where {DT,QobjType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    ψout = QuantumObject(similar(ψ.data), ψ.type, ψ.dims)
+    return A(ψout, ψ, p, t)
+end
+
+@doc raw"""
+    (A::QuantumObjectEvolution)(p, t)
+
+Calculate the time-dependent [`QuantumObjectEvolution`](@ref) object `A` at time `t` with parameters `p`.
+
+# Arguments
+- `p`: The parameters of the time-dependent coefficients.
+- `t`: The time at which the coefficients are evaluated.
+
+# Returns
+- `A::QuantumObject`: The output state.
+"""
+function (A::QuantumObjectEvolution)(p, t)
     # We put 0 in the place of `u` because the time-dependence doesn't depend on the state
-    update_coefficients!(QO.data, 0, p, t)
-    return QuantumObject(concretize(QO.data), QO.type, QO.dims)
+    update_coefficients!(A.data, 0, p, t)
+    return QuantumObject(concretize(A.data), A.type, A.dims)
 end
 
-(QO::QuantumObjectEvolution)(t) = QO((), t)
+(A::QuantumObjectEvolution)(t) = A(nothing, t)
 
-Base.promote_type(A::QuantumObjectEvolution, B::QuantumObjectEvolution) = get_typename_wrapper(A)
-Base.promote_type(A::QuantumObjectEvolution, B::QuantumObject) = get_typename_wrapper(A)
-Base.promote_type(A::QuantumObject, B::QuantumObjectEvolution) = get_typename_wrapper(B)
-Base.promote_type(A::QuantumObject, B::QuantumObject) = get_typename_wrapper(A)
+#=
+`promote_type` should be applied on types. Here I define `promote_op_type` because it is applied to operators.
+=#
+promote_op_type(A::QuantumObjectEvolution, B::QuantumObjectEvolution) = get_typename_wrapper(A)
+promote_op_type(A::QuantumObjectEvolution, B::QuantumObject) = get_typename_wrapper(A)
+promote_op_type(A::QuantumObject, B::QuantumObjectEvolution) = get_typename_wrapper(B)
+promote_op_type(A::QuantumObject, B::QuantumObject) = get_typename_wrapper(A)

src/time_evolution/mcsolve.jl

@@ -81,7 +81,9 @@ function _mcsolve_prob_func(prob, i, repeat)
         ),
     )
 
-    return remake(prob, p = prm)
+    f = deepcopy(prob.f.f)
+
+    return remake(prob, f = f, p = prm)
 end
 
 # Standard output function

src/time_evolution/ssesolve.jl

@@ -175,7 +175,7 @@ function ssesolveProblem(
     check_dims(H_eff_evo, ψ0)
     dims = H_eff_evo.dims
 
-    ψ0 = get_data(ψ0)
+    ψ0 = sparse_to_dense(_CType(ψ0), get_data(ψ0))
 
     progr = ProgressBar(length(tlist), enable = false)
 

test/core-test/quantum_objects.jl

@@ -120,6 +120,24 @@
         @test_throws DimensionMismatch Qobj(ρ_bra.data, type = OperatorBra, dims = 4)
     end
 
+    @testset "Checks on non-QuantumObjects" begin
+        x = 1
+        @test isket(x) == false
+        @test isbra(x) == false
+        @test isoper(x) == false
+        @test issuper(x) == false
+        @test isoperket(x) == false
+        @test isoperbra(x) == false
+
+        x = rand(ComplexF64, 2)
+        @test isket(x) == false
+        @test isbra(x) == false
+        @test isoper(x) == false
+        @test issuper(x) == false
+        @test isoperket(x) == false
+        @test isoperbra(x) == false
+    end
+
     @testset "arithmetic" begin
         a = sprand(ComplexF64, 100, 100, 0.1)
         a2 = Qobj(a)

test/core-test/quantum_objects_evo.jl

@@ -1,17 +1,24 @@
 @testset "Quantum Objects Evolution" verbose = true begin
     # DomainError: incompatible between size of array and type
-    @testset "DomainError" begin
+    @testset "Thrown Errors" begin
         a = MatrixOperator(rand(ComplexF64, 3, 2))
         for t in [Operator, SuperOperator]
             @test_throws DomainError QobjEvo(a, type = t)
         end
-    end
 
-    @testset "ArgumentErro" begin
         a = MatrixOperator(rand(ComplexF64, 3, 2))
         for t in (Ket, Bra, OperatorKet, OperatorBra)
             @test_throws ArgumentError QobjEvo(a, type = t)
         end
+
+        a = QobjEvo(destroy(20))
+        @test_throws ArgumentError QobjEvo(a, type = SuperOperator)
+
+        a = MatrixOperator(rand(ComplexF64, 5, 5))
+        @test_throws DimensionMismatch QobjEvo(a, type = SuperOperator)
+
+        ψ = fock(10, 3)
+        @test_throws TypeError QobjEvo(ψ)
     end
 
     # unsupported type of dims
@@ -47,6 +54,15 @@
         @test_throws DimensionMismatch QobjEvo(a, dims = 2)
     end
 
+    @testset "Promote Operators Type" begin
+        a = destroy(20)
+        A = QobjEvo(a)
+        @test QuantumToolbox.promote_op_type(a, A) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(A, a) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(A, A) == QuantumObjectEvolution
+        @test QuantumToolbox.promote_op_type(a, a) == QuantumObject
+    end
+
     @testset "arithmetic" begin
         a = MatrixOperator(sprand(ComplexF64, 100, 100, 0.1))
         a2 = QobjEvo(a)
@@ -65,7 +81,7 @@
         @test adjoint(a2).data == adjoint(a2.data)
 
         N = 10
-        a = QobjEvo(destroy(N))
+        a = QobjEvo(destroy(N), Operator, N)
         a_d = a'
         X = a + a_d
         # Y = 1im * (a - a_d) # Currently doesn't work. Fix in SciMLOperators.jl
@@ -148,5 +164,16 @@
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])
 
         op1 = QobjEvo(((a, coef1), a' * a, (a', coef2)))
+        op1 = QobjEvo(((a, coef1), a' * a, (a', coef2)))
+
+        p = (ω1 = 1, ω2 = 2)
+        @test op1(p, 0.1) ≈ coef1(p, 0.1) * a + a' * a + coef2(p, 0.1) * a'
+
+        ψ = fock(N, 1)
+        @test op1(ψ, p, 0.1) ≈ (coef1(p, 0.1) * a + a' * a + coef2(p, 0.1) * a') * ψ
+
+        @test iscontant(a) == true
+        @test iscontant(op1) == false
+        @test isconstant(Qobj(a)) == true
     end
 end