Introduce QobjEvo and use SciMLOperators for time evolution (#266)

* First working implementation

* Minor changes

* First working case of sesolve

* Minor changes

* Rebase commits

* Apply Yi-Te comments

* Working mesolve

* Make dsf_mesolve and dfd_mesolve work

* Minor changes

* Working version of ssesolve

* Remove sleep in runtests

* Remove OperatorSum and TimeDependentOperatorSum

* Remove alg as argument from all problem definitions

* First working  runtests

* Add tests for QObjEvo and time evolution

* Add docstrings to documentation

* Reduce tolerance for stochastic runtests

* Make runtests multithreaded and reduce time for timeevolution tests

* Add QobjEvo tests and minor changes

* Update docstrings

* Add comment on the use of MatrixOperator

Author: albertomercurio

.github/workflows/CI.yml

@@ -71,6 +71,7 @@ jobs:
       - uses: julia-actions/julia-runtest@v1
         env:
           GROUP: ${{ matrix.group }}
+          JULIA_NUM_THREADS: auto
       - uses: julia-actions/julia-processcoverage@v1
         with:
           directories: src,ext

docs/src/api.md

@@ -13,6 +13,7 @@ Pages = ["api.md"]
 ## [Quantum object (Qobj) and type](@id doc-API:Quantum-object-and-type)
 
 ```@docs
+AbstractQuantumObject
 BraQuantumObject
 Bra
 KetQuantumObject
@@ -26,10 +27,15 @@ OperatorKet
 SuperOperatorQuantumObject
 SuperOperator
 QuantumObject
-OperatorSum
+QuantumObjectEvolution
 size
 eltype
 length
+```
+
+## [Qobj boolean functions](@id doc-API:Qobj-boolean-functions)
+
+```@docs
 isbra
 isket
 isoper
@@ -40,6 +46,7 @@ LinearAlgebra.ishermitian
 LinearAlgebra.issymmetric
 LinearAlgebra.isposdef
 isunitary
+isconstant
 ```
 
 ## [Qobj arithmetic and attributes](@id doc-API:Qobj-arithmetic-and-attributes)
@@ -154,6 +161,7 @@ lindblad_dissipator
 ## [Synonyms of functions for Qobj](@id doc-API:Synonyms-of-functions-for-Qobj)
 ```@docs
 Qobj
+QobjEvo
 shape
 isherm
 trans

src/QuantumToolbox.jl

@@ -21,6 +21,7 @@ import SciMLBase:
     reinit!,
     remake,
     u_modified!,
+    ODEFunction,
     ODEProblem,
     SDEProblem,
     EnsembleProblem,
@@ -32,13 +33,21 @@ import SciMLBase:
     ContinuousCallback,
     DiscreteCallback
 import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1
-import SciMLOperators: MatrixOperator
+import SciMLOperators:
+    AbstractSciMLOperator,
+    MatrixOperator,
+    ScalarOperator,
+    IdentityOperator,
+    cache_operator,
+    update_coefficients!,
+    concretize,
+    isconstant
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
 import OrdinaryDiffEqTsit5: Tsit5
-import DiffEqNoiseProcess: RealWienerProcess
+import DiffEqNoiseProcess: RealWienerProcess!
 
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
@@ -62,7 +71,9 @@ include("progress_bar.jl")
 include("linear_maps.jl")
 
 # Quantum Object
+include("qobj/quantum_object_base.jl")
 include("qobj/quantum_object.jl")
+include("qobj/quantum_object_evo.jl")
 include("qobj/boolean_functions.jl")
 include("qobj/arithmetic_and_attributes.jl")
 include("qobj/eigsolve.jl")
@@ -71,7 +82,6 @@ include("qobj/states.jl")
 include("qobj/operators.jl")
 include("qobj/superoperators.jl")
 include("qobj/synonyms.jl")
-include("qobj/operator_sum.jl")
 
 # time evolution
 include("time_evolution/time_evolution.jl")

src/qobj/arithmetic_and_attributes.jl

@@ -36,84 +36,80 @@ end
 
 for op in (:(+), :(-), :(*))
     @eval begin
-        function LinearAlgebra.$op(
-            A::QuantumObject{<:AbstractArray{T1},OpType},
-            B::QuantumObject{<:AbstractArray{T2},OpType},
-        ) where {T1,T2,OpType<:QuantumObjectType}
-            A.dims != B.dims &&
-                throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
-            return QuantumObject($(op)(A.data, B.data), A.type, A.dims)
+        function LinearAlgebra.$op(A::AbstractQuantumObject, B::AbstractQuantumObject)
+            check_dims(A, B)
+            QType = promote_op_type(A, B)
+            return QType($(op)(A.data, B.data), A.type, A.dims)
         end
-        LinearAlgebra.$op(A::QuantumObject{<:AbstractArray{T}}) where {T} = QuantumObject($(op)(A.data), A.type, A.dims)
+        LinearAlgebra.$op(A::AbstractQuantumObject) = get_typename_wrapper(A)($(op)(A.data), A.type, A.dims)
 
-        LinearAlgebra.$op(n::T1, A::QuantumObject{<:AbstractArray{T2}}) where {T1<:Number,T2} =
-            QuantumObject($(op)(n * I, A.data), A.type, A.dims)
-        LinearAlgebra.$op(A::QuantumObject{<:AbstractArray{T1}}, n::T2) where {T1,T2<:Number} =
-            QuantumObject($(op)(A.data, n * I), A.type, A.dims)
+        LinearAlgebra.$op(n::T, A::AbstractQuantumObject) where {T<:Number} =
+            get_typename_wrapper(A)($(op)(n * I, A.data), A.type, A.dims)
+        LinearAlgebra.$op(A::AbstractQuantumObject, n::T) where {T<:Number} =
+            get_typename_wrapper(A)($(op)(A.data, n * I), A.type, A.dims)
     end
 end
 
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,OperatorQuantumObject},
+    B::QuantumObject{DT2,KetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Ket, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},BraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,BraQuantumObject},
+    B::AbstractQuantumObject{DT2,OperatorQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Bra, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},KetQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},BraQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,KetQuantumObject},
+    B::QuantumObject{DT2,BraQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, Operator, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},BraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,BraQuantumObject},
+    B::QuantumObject{DT2,KetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
+    B::QuantumObject{DT2,OperatorQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(vec2mat(A.data * mat2vec(B.data)), Operator, A.dims)
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},OperatorBraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorKetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::QuantumObject{DT1,OperatorBraQuantumObject},
+    B::QuantumObject{DT2,OperatorKetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return A.data * B.data
 end
 function LinearAlgebra.:(*)(
-    A::QuantumObject{<:AbstractArray{T1},SuperOperatorQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},OperatorKetQuantumObject},
-) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    A::AbstractQuantumObject{DT1,SuperOperatorQuantumObject},
+    B::QuantumObject{DT2,OperatorKetQuantumObject},
+) where {DT1,DT2}
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, OperatorKet, A.dims)
 end
 function LinearAlgebra.:(*)(
     A::QuantumObject{<:AbstractArray{T1},OperatorBraQuantumObject},
-    B::QuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
+    B::AbstractQuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
 ) where {T1,T2}
-    A.dims != B.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+    check_dims(A, B)
     return QuantumObject(A.data * B.data, OperatorBra, A.dims)
 end
 
-LinearAlgebra.:(^)(A::QuantumObject{<:AbstractArray{T}}, n::T1) where {T,T1<:Number} =
-    QuantumObject(^(A.data, n), A.type, A.dims)
-LinearAlgebra.:(/)(A::QuantumObject{<:AbstractArray{T}}, n::T1) where {T,T1<:Number} =
-    QuantumObject(/(A.data, n), A.type, A.dims)
+LinearAlgebra.:(^)(A::QuantumObject{DT}, n::T) where {DT,T<:Number} = QuantumObject(^(A.data, n), A.type, A.dims)
+LinearAlgebra.:(/)(A::AbstractQuantumObject{DT}, n::T) where {DT,T<:Number} =
+    get_typename_wrapper(A)(A.data / n, A.type, A.dims)
 
 @doc raw"""
     dot(A::QuantumObject, B::QuantumObject)
@@ -125,93 +121,91 @@ Note that `A` and `B` should be [`Ket`](@ref) or [`OperatorKet`](@ref)
 `A ⋅ B` (where `⋅` can be typed by tab-completing `\cdot` in the REPL) is a synonym for `dot(A, B)`
 """
 function LinearAlgebra.dot(
-    A::QuantumObject{<:AbstractArray{T1},OpType},
-    B::QuantumObject{<:AbstractArray{T2},OpType},
-) where {T1<:Number,T2<:Number,OpType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
+    A::QuantumObject{DT1,OpType},
+    B::QuantumObject{DT2,OpType},
+) where {DT1,DT2,OpType<:Union{KetQuantumObject,OperatorKetQuantumObject}}
     A.dims != B.dims && throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(A.data, B.data)
 end
 
 @doc raw"""
-    dot(i::QuantumObject, A::QuantumObject j::QuantumObject)
+    dot(i::QuantumObject, A::AbstractQuantumObject j::QuantumObject)
 
-Compute the generalized dot product `dot(i, A*j)` between three [`QuantumObject`](@ref): ``\langle i | \hat{A} | j \rangle``
+Compute the generalized dot product `dot(i, A*j)` between a [`AbstractQuantumObject`](@ref) and two [`QuantumObject`](@ref) (`i` and `j`), namely ``\langle i | \hat{A} | j \rangle``.
 
 Supports the following inputs:
 - `A` is in the type of [`Operator`](@ref), with `i` and `j` are both [`Ket`](@ref).
 - `A` is in the type of [`SuperOperator`](@ref), with `i` and `j` are both [`OperatorKet`](@ref)
 """
 function LinearAlgebra.dot(
-    i::QuantumObject{<:AbstractArray{T1},KetQuantumObject},
-    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
-    j::QuantumObject{<:AbstractArray{T3},KetQuantumObject},
-) where {T1<:Number,T2<:Number,T3<:Number}
+    i::QuantumObject{DT1,KetQuantumObject},
+    A::AbstractQuantumObject{DT2,OperatorQuantumObject},
+    j::QuantumObject{DT3,KetQuantumObject},
+) where {DT1,DT2,DT3}
     ((i.dims != A.dims) || (A.dims != j.dims)) &&
         throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 function LinearAlgebra.dot(
-    i::QuantumObject{<:AbstractArray{T1},OperatorKetQuantumObject},
-    A::QuantumObject{<:AbstractArray{T2},SuperOperatorQuantumObject},
-    j::QuantumObject{<:AbstractArray{T3},OperatorKetQuantumObject},
-) where {T1<:Number,T2<:Number,T3<:Number}
+    i::QuantumObject{DT1,OperatorKetQuantumObject},
+    A::AbstractQuantumObject{DT2,SuperOperatorQuantumObject},
+    j::QuantumObject{DT3,OperatorKetQuantumObject},
+) where {DT1,DT2,DT3}
     ((i.dims != A.dims) || (A.dims != j.dims)) &&
         throw(DimensionMismatch("The quantum objects are not of the same Hilbert dimension."))
     return LinearAlgebra.dot(i.data, A.data, j.data)
 end
 
 @doc raw"""
-    conj(A::QuantumObject)
+    conj(A::AbstractQuantumObject)
 
-Return the element-wise complex conjugation of the [`QuantumObject`](@ref).
+Return the element-wise complex conjugation of the [`AbstractQuantumObject`](@ref).
 """
-Base.conj(A::QuantumObject{<:AbstractArray{T}}) where {T} = QuantumObject(conj(A.data), A.type, A.dims)
+Base.conj(A::AbstractQuantumObject) = get_typename_wrapper(A)(conj(A.data), A.type, A.dims)
 
 @doc raw"""
-    transpose(A::QuantumObject)
+    transpose(A::AbstractQuantumObject)
 
-Lazy matrix transpose of the [`QuantumObject`](@ref).
+Lazy matrix transpose of the [`AbstractQuantumObject`](@ref).
 """
 LinearAlgebra.transpose(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(transpose(A.data), A.type, A.dims)
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    get_typename_wrapper(A)(transpose(A.data), A.type, A.dims)
 
 @doc raw"""
     A'
-    adjoint(A::QuantumObject)
+    adjoint(A::AbstractQuantumObject)
 
-Lazy adjoint (conjugate transposition) of the [`QuantumObject`](@ref)
+Lazy adjoint (conjugate transposition) of the [`AbstractQuantumObject`](@ref)
 
 Note that `A'` is a synonym for `adjoint(A)`
 """
 LinearAlgebra.adjoint(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
-    QuantumObject(adjoint(A.data), A.type, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},KetQuantumObject}) where {T} =
-    QuantumObject(adjoint(A.data), Bra, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},BraQuantumObject}) where {T} =
-    QuantumObject(adjoint(A.data), Ket, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},OperatorKetQuantumObject}) where {T} =
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    get_typename_wrapper(A)(adjoint(A.data), A.type, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,KetQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Bra, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,BraQuantumObject}) where {DT} = QuantumObject(adjoint(A.data), Ket, A.dims)
+LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorKetQuantumObject}) where {DT} =
     QuantumObject(adjoint(A.data), OperatorBra, A.dims)
-LinearAlgebra.adjoint(A::QuantumObject{<:AbstractArray{T},OperatorBraQuantumObject}) where {T} =
+LinearAlgebra.adjoint(A::QuantumObject{DT,OperatorBraQuantumObject}) where {DT} =
     QuantumObject(adjoint(A.data), OperatorKet, A.dims)
 
 @doc raw"""
-    inv(A::QuantumObject)
+    inv(A::AbstractQuantumObject)
 
-Matrix inverse of the [`QuantumObject`](@ref)
+Matrix inverse of the [`AbstractQuantumObject`](@ref). If `A` is a [`QuantumObjectEvolution`](@ref), the inverse is computed at the last computed time.
 """
 LinearAlgebra.inv(
-    A::QuantumObject{<:AbstractArray{T},OpType},
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+    A::AbstractQuantumObject{DT,OpType},
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     QuantumObject(sparse(inv(Matrix(A.data))), A.type, A.dims)
 
 LinearAlgebra.Hermitian(
-    A::QuantumObject{<:AbstractArray{T},OpType},
+    A::QuantumObject{DT,OpType},
     uplo::Symbol = :U,
-) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
+) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} =
     QuantumObject(Hermitian(A.data, uplo), A.type, A.dims)
 
 @doc raw"""
@@ -436,8 +430,8 @@ Matrix sine of [`QuantumObject`](@ref), defined as
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
 LinearAlgebra.sin(
-    A::QuantumObject{<:AbstractMatrix{T},ObjType},
-) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) - exp(-1im * A)) / 2im
+    A::QuantumObject{DT,ObjType},
+) where {DT,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) - exp(-1im * A)) / 2im
 
 @doc raw"""
     cos(A::QuantumObject)
@@ -449,8 +443,8 @@ Matrix cosine of [`QuantumObject`](@ref), defined as
 Note that this function only supports for [`Operator`](@ref) and [`SuperOperator`](@ref)
 """
 LinearAlgebra.cos(
-    A::QuantumObject{<:AbstractMatrix{T},ObjType},
-) where {T,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) + exp(-1im * A)) / 2
+    A::QuantumObject{DT,ObjType},
+) where {DT,ObjType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = (exp(1im * A) + exp(-1im * A)) / 2
 
 @doc raw"""
     diag(A::QuantumObject, k::Int=0)
@@ -659,11 +653,11 @@ tidyup!(A::AbstractArray{T}, tol::T2 = 1e-14) where {T,T2<:Real} =
     @. A = T(abs(real(A)) > tol) * real(A) + 1im * T(abs(imag(A)) > tol) * imag(A)
 
 @doc raw"""
-    get_data(A::QuantumObject)
+    get_data(A::AbstractQuantumObject)
 
-Returns the data of a QuantumObject.
+Returns the data of a [`AbstractQuantumObject`](@ref).
 """
-get_data(A::QuantumObject) = A.data
+get_data(A::AbstractQuantumObject) = A.data
 
 @doc raw"""
     get_coherence(ψ::QuantumObject)