Working mesolve

Author: albertomercurio

src/qobj/boolean_functions.jl

@@ -93,6 +93,6 @@ isunitary(U::QuantumObject{<:AbstractArray{T}}; kwargs...) where {T} =
 @doc raw"""
     isconstant(A::AbstractQuantumObject)
 
-Test whether the [`AbstractQuantumObject`](@ref) `A` is constant in time. For a [`QuantumObject`](@ref), this function returns `false`, while for a [`QuantumObjectEvolution`](@ref), this function returns `true` if the operator is contant in time.
+Test whether the [`AbstractQuantumObject`](@ref) `A` is constant in time. For a [`QuantumObject`](@ref), this function returns `true`, while for a [`QuantumObjectEvolution`](@ref), this function returns `true` if the operator is contant in time.
 """
 isconstant(A::AbstractQuantumObject) = isconstant(A.data)

src/qobj/eigsolve.jl

@@ -362,26 +362,24 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 - [1] Minganti, F., & Huybrechts, D. (2022). Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of time-independent and Floquet open quantum systems. Quantum, 6, 649.
 """
 function eigsolve_al(
-    H::QuantumObject{MT1,HOpType},
+    H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
     T::Real,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
-    H_t::Union{Nothing,Function} = nothing,
     params::NamedTuple = NamedTuple(),
     ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
     k::Int = 1,
     krylovdim::Int = min(10, size(H, 1)),
     maxiter::Int = 200,
     eigstol::Real = 1e-6,
     kwargs...,
-) where {MT1<:AbstractMatrix,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
-    L = liouvillian(H, c_ops)
+) where {DT1,HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    L_evo = _mesolve_make_L_QobjEvo(H, c_ops)
     prob = mesolveProblem(
-        L,
+        L_evo,
         QuantumObject(ρ0, type = Operator, dims = H.dims),
         [zero(T), T];
         alg = alg,
-        H_t = H_t,
         params = params,
         progress_bar = Val(false),
         kwargs...,
@@ -390,17 +388,17 @@ function eigsolve_al(
 
     # prog = ProgressUnknown(desc="Applications:", showspeed = true, enabled=progress)
 
-    Lmap = ArnoldiLindbladIntegratorMap(eltype(MT1), size(L), integrator)
+    Lmap = ArnoldiLindbladIntegratorMap(eltype(DT1), size(L_evo), integrator)
 
-    res = _eigsolve(Lmap, mat2vec(ρ0), L.type, L.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
+    res = _eigsolve(Lmap, mat2vec(ρ0), L_evo.type, L_evo.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
     # finish!(prog)
 
     vals = similar(res.values)
     vecs = similar(res.vectors)
 
     for i in eachindex(res.values)
         vec = view(res.vectors, :, i)
-        vals[i] = dot(vec, L.data, vec)
+        vals[i] = dot(vec, L_evo.data, vec)
         @. vecs[:, i] = vec * exp(-1im * angle(vec[1]))
     end
 

src/qobj/quantum_object.jl

@@ -1,9 +1,6 @@
 #=
-This file defines:
-    1. the QuantumObject (Qobj) structure
-    2. all the type structures for QuantumObject
-Also support for fundamental functions in Julia standard library:
-    - Base: show, length, size, eltype, getindex, setindex!, isequal, :(==), isapprox, Vector, Matrix
+This file defines the QuantumObject (Qobj) structure.
+It also implements the fundamental functions in Julia standard library:
     - SparseArrays: sparse, nnz, nonzeros, rowvals, droptol!, dropzeros, dropzeros!, SparseVector, SparseMatrixCSC
 =#
 

src/qobj/quantum_object_base.jl

@@ -1,3 +1,9 @@
+#=
+This file defines the AbstractQuantumObject structure, all the type structures for QuantumObject and fundamental functions in Julia standard library.
+Also support for fundamental functions in Julia standard library:
+    - Base: show, length, size, eltype, getindex, setindex!, isequal, :(==), isapprox, Vector, Matrix
+=#
+
 export AbstractQuantumObject
 export QuantumObjectType,
     BraQuantumObject,

src/qobj/quantum_object_evo.jl

@@ -82,6 +82,7 @@ Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=fal
 ⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
 ⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
 ⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+````
 """
 struct QuantumObjectEvolution{
     DT<:AbstractSciMLOperator,
@@ -94,8 +95,8 @@ struct QuantumObjectEvolution{
 end
 
 # Make the QuantumObjectEvolution, with the option to pre-multiply by a scalar
-function QuantumObjectEvolution(op_func_list::Tuple, α::Union{Nothing,Number} = nothing)
-    op, data = _generate_data(op_func_list, α)
+function QuantumObjectEvolution(op_func_list::Tuple, α::Union{Nothing,Number} = nothing; f::Function = identity)
+    op, data = _QobjEvo_generate_data(op_func_list, α; f = f)
     dims = op.dims
     type = op.type
 
@@ -106,23 +107,35 @@ function QuantumObjectEvolution(op_func_list::Tuple, α::Union{Nothing,Number} =
     return QuantumObjectEvolution(data, type, dims)
 end
 
-QuantumObjectEvolution(op::QuantumObject, α::Union{Nothing,Number} = nothing) =
+QuantumObjectEvolution(op::QuantumObject, α::Union{Nothing,Number} = nothing; f::Function = identity) =
     QuantumObjectEvolution(_make_SciMLOperator(op, α), op.type, op.dims)
 
-function QuantumObjectEvolution(op::QuantumObjectEvolution, α::Union{Nothing,Number} = nothing)
+function QuantumObjectEvolution(op::QuantumObjectEvolution, α::Union{Nothing,Number} = nothing; f::Function = identity)
+    f !== identity && throw(ArgumentError("The function `f` is not supported for QuantumObjectEvolution inputs."))
     if α isa Nothing
         return QuantumObjectEvolution(op.data, op.type, op.dims)
     end
     return QuantumObjectEvolution(α * op.data, op.type, op.dims)
 end
 
-@generated function _generate_data(op_func_list::Tuple, α)
+#=
+    _QobjEvo_generate_data(op_func_list::Tuple, α; f::Function=identity)
+
+Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution` object. The `op_func_list` is a tuple of tuples or operators. Each element of the tuple can be a tuple with two elements (operator, function) or an operator. The function is used to generate the time-dependent coefficients for the operators. The `α` parameter is used to pre-multiply the operators by a scalar. The `f` parameter is used to pre-applying a function to the operators before converting them to SciML operators. During the parsing, the dimensions of the operators are checked to be the same, and all the constant operators are summed together.
+
+# Arguments
+- `op_func_list::Tuple`: A tuple of tuples or operators.
+- `α`: A scalar to pre-multiply the operators.
+- `f::Function=identity`: A function to pre-apply to the operators.
+=#
+@generated function _QobjEvo_generate_data(op_func_list::Tuple, α; f::Function = identity)
     op_func_list_types = op_func_list.parameters
     N = length(op_func_list_types)
 
     dims_expr = ()
     first_op = nothing
     data_expr = :(0)
+    qobj_expr_const = :(0)
 
     for i in 1:N
         op_func_type = op_func_list_types[i]
@@ -136,49 +149,54 @@ end
                 ),
             )
 
+            op = :(op_func_list[$i][1])
             data_type = op_type.parameters[1]
-            dims_expr = (dims_expr..., :(op_func_list[$i][1].dims))
+            dims_expr = (dims_expr..., :($op.dims))
             if i == 1
-                first_op = :(op_func_list[$i][1])
+                first_op = :(f($op))
             end
+            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α, f = f))
         else
             op_type = op_func_type
             (isoper(op_type) || issuper(op_type)) ||
                 throw(ArgumentError("The element must be a Operator or SuperOperator."))
 
             data_type = op_type.parameters[1]
             dims_expr = (dims_expr..., :(op_func_list[$i].dims))
-
             if i == 1
                 first_op = :(op_func_list[$i])
             end
+            qobj_expr_const = :($qobj_expr_const + f(op_func_list[$i]))
         end
-        data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α))
     end
 
+    data_expr_const = :(_make_SciMLOperator($qobj_expr_const, α))
+
     quote
         dims = tuple($(dims_expr...))
 
         length(unique(dims)) == 1 || throw(ArgumentError("The dimensions of the operators must be the same."))
 
-        return $first_op, $data_expr
+        data_expr = $data_expr_const + $data_expr
+
+        return $first_op, data_expr
     end
 end
 
-function _make_SciMLOperator(op_func::Tuple, α)
+function _make_SciMLOperator(op_func::Tuple, α; f::Function = identity)
     T = eltype(op_func[1])
     update_func = (a, u, p, t) -> op_func[2](p, t)
     if α isa Nothing
-        return ScalarOperator(zero(T), update_func) * MatrixOperator(op_func[1].data)
+        return ScalarOperator(zero(T), update_func) * MatrixOperator(f(op_func[1]).data)
     end
-    return ScalarOperator(zero(T), update_func) * MatrixOperator(α * op_func[1].data)
+    return ScalarOperator(zero(T), update_func) * MatrixOperator(α * f(op_func[1]).data)
 end
 
-function _make_SciMLOperator(op::QuantumObject, α)
+function _make_SciMLOperator(op::QuantumObject, α; f::Function = identity)
     if α isa Nothing
-        return MatrixOperator(op.data)
+        return MatrixOperator(f(op).data)
     end
-    return MatrixOperator(α * op.data)
+    return MatrixOperator(α * f(op.data))
 end
 
 function (QO::QuantumObjectEvolution)(p, t)

src/qobj/synonyms.jl

@@ -18,14 +18,91 @@ Note that this functions is same as `QuantumObject(A; type=type, dims=dims)`.
 Qobj(A; kwargs...) = QuantumObject(A; kwargs...)
 
 @doc raw"""
-    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing)
+    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; f::Function=identity)
 
 Generate [`QuantumObjectEvolution`](@ref)
 
-Note that this functions is same as `QuantumObjectEvolution(op_func_list)`. If `α` is provided, all the operators in `op_func_list` will be pre-multiplied by `α`.
+Note that this functions is same as `QuantumObjectEvolution(op_func_list)`. If `α` is provided, all the operators in `op_func_list` will be pre-multiplied by `α`. The `f` parameter is used to pre-apply a function to the operators before converting them to SciML operators.
+
+# Arguments
+- `op_func_list::Union{Tuple,AbstractQuantumObject}`: A tuple of tuples or operators.
+- `α::Union{Nothing,Number}=nothing`: A scalar to pre-multiply the operators.
+- `f::Function=identity`: A function to pre-apply to the operators.
+
+!!! warning "Beware of type-stability!"
+    Please note that, unlike QuTiP, this function doesn't support `op_func_list` as `Vector` type. This is related to the type-stability issue. See the Section [The Importance of Type-Stability](@ref doc:Type-Stability) for more details.
+
+# Examples
+This operator can be initialized in the same way as the QuTiP `QobjEvo` object. For example
+```
+julia> a = tensor(destroy(10), qeye(2))
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 18 stored entries:
+⎡⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦
+
+julia> σm = tensor(qeye(10), sigmam())
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
+⎡⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦
+
+julia> coef1(p, t) = exp(-1im * t)
+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)))
+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))
+```
+
+We can also concretize the operator at a specific time `t`
+```
+julia> op1(0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+```
+
+It also supports parameter-dependent time evolution
+```
+julia> coef1(p, t) = exp(-1im * p.ω1 * t)
+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)))
+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))
+
+julia> p = (ω1 = 1.0, ω2 = 0.5)
+(ω1 = 1.0, ω2 = 0.5)
+
+julia> op1(p, 0.1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
+20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
+⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
+⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
+⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
+````
 """
-QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number} = nothing) =
-    QuantumObjectEvolution(op_func_list, α)
+QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number} = nothing; f::Function = identity) =
+    QuantumObjectEvolution(op_func_list, α; f = f)
 
 @doc raw"""
     shape(A::QuantumObject)

src/steadystate.jl

@@ -179,15 +179,15 @@ _steadystate(
     kwargs...,
 ) where {T} = throw(
     ArgumentError(
-        "The initial state ψ0 is required for SteadyStateODESolver, use the following call instead: `steadystate(H, ψ0, tspan, c_ops)`.",
+        "The initial state ψ0 is required for SteadyStateODESolver, use the following call instead: `steadystate(H, ψ0, tmax, c_ops)`.",
     ),
 )
 
 @doc raw"""
     steadystate(
         H::QuantumObject,
         ψ0::QuantumObject,
-        tspan::Real = Inf,
+        tmax::Real = Inf,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateODESolver = SteadyStateODESolver(),
         reltol::Real = 1.0e-8,
@@ -212,48 +212,34 @@ or
 # Parameters
 - `H::QuantumObject`: The Hamiltonian or the Liouvillian of the system.
 - `ψ0::QuantumObject`: The initial state of the system.
-- `tspan::Real=Inf`: The final time step for the steady state problem.
+- `tmax::Real=Inf`: The final time step for the steady state problem.
 - `c_ops::Union{Nothing,AbstractVector,Tuple}=nothing`: The list of the collapse operators.
 - `solver::SteadyStateODESolver=SteadyStateODESolver()`: see [`SteadyStateODESolver`](@ref) for more details.
 - `reltol::Real=1.0e-8`: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
 - `abstol::Real=1.0e-10`: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
 - `kwargs...`: The keyword arguments for the ODEProblem.
 """
 function steadystate(
-    H::QuantumObject{MT1,HOpType},
-    ψ0::QuantumObject{<:AbstractArray{T2},StateOpType},
-    tspan::Real = Inf,
+    H::QuantumObject{DT1,HOpType},
+    ψ0::QuantumObject{DT2,StateOpType},
+    tmax::Real = Inf,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateODESolver = SteadyStateODESolver(),
     reltol::Real = 1.0e-8,
     abstol::Real = 1.0e-10,
     kwargs...,
 ) where {
-    MT1<:AbstractMatrix,
-    T2,
+    DT1,
+    DT2,
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject},
 }
-    check_dims(H, ψ0)
-
-    N = prod(H.dims)
-    u0 = sparse_to_dense(_CType(ψ0), mat2vec(ket2dm(ψ0).data))
-
-    L = MatrixOperator(liouvillian(H, c_ops).data)
-
     ftype = _FType(ψ0)
-    prob = ODEProblem{true}(L, u0, (ftype(0), ftype(tspan))) # Convert tspan to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
-    sol = solve(
-        prob,
-        solver.alg;
-        callback = TerminateSteadyState(abstol, reltol, _steadystate_ode_condition),
-        reltol = reltol,
-        abstol = abstol,
-        kwargs...,
-    )
+    cb = TerminateSteadyState(abstol, reltol, _steadystate_ode_condition)
+    sol = mesolve(H, ψ0, [ftype(0), ftype(tmax)], c_ops, progress_bar = Val(false), callback = cb)
 
-    ρss = reshape(sol.u[end], N, N)
-    return QuantumObject(ρss, Operator, H.dims)
+    ρss = sol.states[end]
+    return ρss
 end
 
 function _steadystate_ode_condition(integrator, abstol, reltol, min_t)