A lazy method to construct QobjEvo (qutip#282)

* add a lazy method to construct `QobjEvo`

* format files

* improve tests to make sure time-independent problems are `MatrixOperator`

* modify docstrings

* fix typo

* fix extra `ScalarOperator` in `mcsolve`

* rebase to `ScaledOperator` in time-independent `mcsolve`

* add missing import in test

Author: ytdHuang

src/qobj/quantum_object_evo.jl

@@ -27,6 +27,17 @@ Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=fal
 ⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
 ⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦
 
+julia> coef1(p, t) = exp(-1im * t)
+coef1 (generic function with 1 method)
+
+julia> op = QuantumObjectEvolution(a, coef1)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20)
+```
+
+If there are more than 2 operators, we need to put each set of operator and coefficient function into a two-element `Tuple`, and put all these `Tuple`s together in a larger `Tuple`:
+
+```
 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:
@@ -36,9 +47,6 @@ Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=fal
 ⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
 ⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦
 
-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)
 
@@ -134,7 +142,7 @@ function QuantumObjectEvolution(data::AbstractSciMLOperator, type::QuantumObject
     return QuantumObjectEvolution(data, type, SVector{1,Int}(dims))
 end
 
-"""
+@doc raw"""
     QuantumObjectEvolution(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)
 
 Generate a [`QuantumObjectEvolution`](@ref) object from a [`SciMLOperator`](https://github.com/SciML/SciMLOperators.jl), in the same way as [`QuantumObject`](@ref) for `AbstractArray` inputs.
@@ -172,6 +180,9 @@ function QuantumObjectEvolution(
     return QuantumObjectEvolution(data, type, dims)
 end
 
+QuantumObjectEvolution(op::QuantumObject, f::Function; type::Union{Nothing,QuantumObjectType} = nothing) =
+    QuantumObjectEvolution(((op, f),); type = type)
+
 function QuantumObjectEvolution(
     op::QuantumObject,
     α::Union{Nothing,Number} = nothing;

src/qobj/synonyms.jl

@@ -17,6 +17,35 @@ Note that this functions is same as `QuantumObject(A; type=type, dims=dims)`.
 """
 Qobj(A; kwargs...) = QuantumObject(A; kwargs...)
 
+@doc raw"""
+    QobjEvo(op::QuantumObject, f::Function; type::Union{Nothing,QuantumObjectType} = nothing)
+
+Generate [`QuantumObjectEvolution`](@ref).
+
+Note that this functions is same as `QuantumObjectEvolution(op, f; type = type)`. The `f` parameter is used to pre-apply a function to the operators before converting them to SciML operators. The `type` parameter is used to specify the type of the [`QuantumObject`](@ref), either `Operator` or `SuperOperator`.
+
+# Examples
+```
+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> coef(p, t) = exp(-1im * t)
+coef1 (generic function with 1 method)
+
+julia> op = QobjEvo(a, coef)
+Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
+ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20)
+```
+"""
+QobjEvo(op::QuantumObject, f::Function; type::Union{Nothing,QuantumObjectType} = nothing) =
+    QuantumObjectEvolution(op, f; type = type)
+
 @doc raw"""
     QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type::Union{Nothing, QuantumObjectType}=nothing)
 

test/core-test/quantum_objects_evo.jl

@@ -185,7 +185,7 @@
 
         # SuperOperator
         X = a * a'
-        c_op1 = QobjEvo(((a', coef1),))
+        c_op1 = QobjEvo(a', coef1)
         c_op2 = QobjEvo(((a, coef2), (X, coef3)))
         c_ops = [c_op1, c_op2]
         D1_ti = abs2(coef1(p, t)) * lindblad_dissipator(a')
@@ -213,7 +213,7 @@
         @test_throws ArgumentError cache_operator(L_td, ψ)
 
         @testset "Type Inference" begin
-            @inferred liouvillian(H_td, (a, QobjEvo(((a', coef1),))))
+            @inferred liouvillian(H_td, (a, QobjEvo(a', coef1)))
         end
     end
 end

test/core-test/time_evolution.jl

@@ -11,10 +11,12 @@
         psi0 = kron(fock(N, 0), fock(2, 0))
         t_l = LinRange(0, 1000, 1000)
         e_ops = [a_d * a]
-        sol = sesolve(H, psi0, t_l, e_ops = e_ops, progress_bar = Val(false))
+        prob = sesolveProblem(H, psi0, t_l, e_ops = e_ops, progress_bar = Val(false))
+        sol = sesolve(prob)
         sol2 = sesolve(H, psi0, t_l, progress_bar = Val(false))
         sol3 = sesolve(H, psi0, t_l, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
         sol_string = sprint((t, s) -> show(t, "text/plain", s), sol)
+        @test prob.f.f isa MatrixOperator
         @test sum(abs.(sol.expect[1, :] .- sin.(η * t_l) .^ 2)) / length(t_l) < 0.1
         @test length(sol.times) == length(t_l)
         @test length(sol.states) == 1
@@ -55,9 +57,11 @@
         e_ops = [a_d * a]
         psi0 = basis(N, 3)
         t_l = LinRange(0, 100, 1000)
-        sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
+        prob_me = mesolveProblem(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
+        sol_me = mesolve(prob_me)
         sol_me2 = mesolve(H, psi0, t_l, c_ops, progress_bar = Val(false))
         sol_me3 = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
+        prob_mc = mcsolveProblem(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
         sol_mc = mcsolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false))
         sol_mc2 = mcsolve(
             H,
@@ -93,6 +97,8 @@
         sol_me_string = sprint((t, s) -> show(t, "text/plain", s), sol_me)
         sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
         sol_sse_string = sprint((t, s) -> show(t, "text/plain", s), sol_sse)
+        @test prob_me.f.f isa MatrixOperator
+        @test prob_mc.f.f isa ScaledOperator # TODO: can be optimized as MatrixOperator
         @test sum(abs.(sol_mc.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test sum(abs.(sol_mc2.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test sum(abs.(vec(expect_mc_states_mean) .- vec(sol_me.expect))) / length(t_l) < 0.1
@@ -230,7 +236,7 @@
 
         @testset "Type Inference mesolve" begin
             coef(p, t) = exp(-t)
-            ad_t = QobjEvo(((a', coef),))
+            ad_t = QobjEvo(a', coef)
             @inferred mesolveProblem(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolveProblem(H, ψ0, [0, 10], c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolveProblem(

test/runtests.jl

@@ -4,6 +4,7 @@ using QuantumToolbox
 using QuantumToolbox: position, momentum
 using Random
 using SciMLOperators
+import SciMLOperators: ScaledOperator
 
 const GROUP = get(ENV, "GROUP", "All")