From 3dfd9a0e2b29da0f6a3d223a33092f2271ba0bc1 Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Sat, 27 Sep 2025 17:20:12 +0200
Subject: [PATCH 1/5] Remove alpha parameter from QobjEvo definition

---
 src/qobj/quantum_object_evo.jl | 69 ++++++++++++----------------------
 src/time_evolution/mcsolve.jl  |  4 +-
 src/time_evolution/sesolve.jl  |  8 +---
 src/time_evolution/ssesolve.jl |  2 +-
 4 files changed, 27 insertions(+), 56 deletions(-)

diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index b7a165b4f..a7fd3f551 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -175,19 +175,18 @@ function QuantumObjectEvolution(data::AbstractSciMLOperator; type = Operator(),
 end
 
 @doc raw"""
-    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type=nothing)
-    QuantumObjectEvolution(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type=nothing)
+    QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}; type=nothing)
+    QuantumObjectEvolution(op_func_list::Union{Tuple,AbstractQuantumObject}; type=nothing)
 
 Generate [`QuantumObjectEvolution`](@ref).
 
 # Arguments
 - `op_func_list::Union{Tuple,AbstractQuantumObject}`: A tuple of tuples or operators.
-- `α::Union{Nothing,Number}=nothing`: A scalar to pre-multiply 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.
 
-Note that if `α` is provided, all the operators in `op_func_list` will be pre-multiplied by `α`. The `type` parameter is used to specify the type of the [`QuantumObject`](@ref), either `Operator` or `SuperOperator`. 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`.
 
 !!! note
     `QobjEvo` is a synonym of `QuantumObjectEvolution`.
@@ -267,8 +266,8 @@ Quantum Object:   type=Operator()   dims=[10, 2]   size=(20, 20)   ishermitian=f
 ⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦
 ```
 """
-function QuantumObjectEvolution(op_func_list::Tuple, α::Union{Nothing,Number} = nothing; type = nothing)
-    op, data = _QobjEvo_generate_data(op_func_list, α)
+function QuantumObjectEvolution(op_func_list::Tuple; type = nothing)
+    op, data = _QobjEvo_generate_data(op_func_list)
     dims = op.dimensions
     _check_type(type)
 
@@ -284,17 +283,17 @@ function QuantumObjectEvolution(op_func_list::Tuple, α::Union{Nothing,Number} =
 end
 
 # this is a extra method if user accidentally specify `QuantumObjectEvolution( (op, func) )` or `QuantumObjectEvolution( ((op, func)) )`
-QuantumObjectEvolution(op_func::Tuple{QuantumObject,Function}, α::Union{Nothing,Number} = nothing; type = nothing) =
-    QuantumObjectEvolution((op_func,), α; type = type)
+QuantumObjectEvolution(op_func::Tuple{<:QuantumObject,<:Function}; type = nothing) =
+    QuantumObjectEvolution((op_func,); type = type)
 
 @doc raw"""
-    QuantumObjectEvolution(op::QuantumObject, f::Function, α::Union{Nothing,Number}=nothing; type = nothing)
-    QobjEvo(op::QuantumObject, f::Function, α::Union{Nothing,Number}=nothing; type = nothing)
+    QuantumObjectEvolution(op::QuantumObject, f::Function; type = nothing)
+    QobjEvo(op::QuantumObject, f::Function; type = nothing)
 
 Generate [`QuantumObjectEvolution`](@ref).
 
 # Notes
-- 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`.
+- The `f` parameter is time-dependent coefficient that multiplies the operator. The `type` parameter is used to specify the type of the [`QuantumObject`](@ref), either `Operator` or `SuperOperator`.
 - `QobjEvo` is a synonym of `QuantumObjectEvolution`.
 
 # Examples
@@ -318,18 +317,17 @@ Quantum Object Evo.:   type=Operator()   dims=[10, 2]   size=(20, 20)   ishermit
 ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20)
 ```
 """
-QuantumObjectEvolution(op::QuantumObject, f::Function, α::Union{Nothing,Number} = nothing; type = nothing) =
-    QuantumObjectEvolution(((op, f),), α; type = type)
+QuantumObjectEvolution(op::QuantumObject, f::Function; type = nothing) = QuantumObjectEvolution(((op, f),); type = type)
 
-function QuantumObjectEvolution(op::QuantumObject, α::Union{Nothing,Number} = nothing; type = nothing)
+function QuantumObjectEvolution(op::QuantumObject; type = nothing)
     _check_type(type)
     if type isa Nothing
         type = op.type
     end
-    return QuantumObjectEvolution(_make_SciMLOperator(op, α), type, op.dimensions)
+    return QuantumObjectEvolution(_make_SciMLOperator(op), type, op.dimensions)
 end
 
-function QuantumObjectEvolution(op::QuantumObjectEvolution, α::Union{Nothing,Number} = nothing; type = nothing)
+function QuantumObjectEvolution(op::QuantumObjectEvolution; type = nothing)
     _check_type(type)
     if type isa Nothing
         type = op.type
@@ -340,23 +338,18 @@ function QuantumObjectEvolution(op::QuantumObjectEvolution, α::Union{Nothing,Nu
             ),
         )
     end
-    if α isa Nothing
-        return QuantumObjectEvolution(op.data, type, op.dimensions)
-    end
-    return QuantumObjectEvolution(_promote_to_scimloperator(α, op.data), type, op.dimensions)
+    return QuantumObjectEvolution(op.data, type, op.dimensions)
 end
 
 #=
-    _QobjEvo_generate_data(op_func_list::Tuple, α; f::Function=identity)
+    _QobjEvo_generate_data(op_func_list::Tuple)
 
-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.
+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. 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, α)
+@generated function _QobjEvo_generate_data(op_func_list::Tuple)
     op_func_list_types = op_func_list.parameters
     N = length(op_func_list_types)
 
@@ -385,7 +378,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
             if i == 1
                 first_op = :($op)
             end
-            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i], α))
+            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i]))
         else
             op_type = op_func_type
             (isoper(op_type) || issuper(op_type)) ||
@@ -414,7 +407,7 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
             )
         end
 
-        data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const, α)
+        data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const)
 
         data_expr = data_expr_const + $data_expr
 
@@ -422,32 +415,16 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
     end
 end
 
-function _make_SciMLOperator(op_func::Tuple, α)
+function _make_SciMLOperator(op_func::Tuple)
     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) * _promote_to_scimloperator(op_func[1].data)
-    end
-    return ScalarOperator(zero(T), update_func) * _promote_to_scimloperator(α, op_func[1].data)
+    return ScalarOperator(zero(T), update_func) * _promote_to_scimloperator(op_func[1].data)
 end
 
-function _make_SciMLOperator(op::AbstractQuantumObject, α)
-    if α isa Nothing
-        return _promote_to_scimloperator(op.data)
-    end
-    return _promote_to_scimloperator(α, op.data)
-end
+_make_SciMLOperator(op::AbstractQuantumObject) = _promote_to_scimloperator(op.data)
 
 _promote_to_scimloperator(data::AbstractMatrix) = MatrixOperator(data)
 _promote_to_scimloperator(data::AbstractSciMLOperator) = data
-_promote_to_scimloperator(α::Number, data::AbstractMatrix) = MatrixOperator(α * data)
-# We still have to define this for AddedOperator, as it is not present in SciMLOperators.jl
-function _promote_to_scimloperator(α::Number, data::AddedOperator)
-    return AddedOperator(_promote_to_scimloperator.(α, data.ops)) # Try to propagate the rule
-end
-function _promote_to_scimloperator(α::Number, data::AbstractSciMLOperator)
-    return α * data # Going back to the generic case
-end
 
 @doc raw"""
     (A::QuantumObjectEvolution)(ψout, ψin, p, t)
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 396fecafd..98b60076e 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -31,11 +31,11 @@ function _mcsolve_make_Heff_QobjEvo(H::QuantumObject, c_ops)
 end
 function _mcsolve_make_Heff_QobjEvo(H::Tuple, c_ops)
     c_ops isa Nothing && return QobjEvo(H)
-    return QobjEvo((H..., -1im * mapreduce(op -> op' * op, +, c_ops) / 2))
+    return QobjEvo((H..., -1im * sum(op -> op' * op, c_ops) / 2))
 end
 function _mcsolve_make_Heff_QobjEvo(H::QuantumObjectEvolution, c_ops)
     c_ops isa Nothing && return H
-    return H + QobjEvo(mapreduce(op -> op' * op, +, c_ops), -1im / 2)
+    return H + QobjEvo(sum(op -> -1im * op' * op / 2, c_ops))
 end
 
 @doc raw"""
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index f7fd3f7c4..3ad229402 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -1,12 +1,6 @@
 export sesolveProblem, sesolve
 
-_sesolve_make_U_QobjEvo(
-    H::QuantumObjectEvolution{Operator,DimsType,<:MatrixOperator},
-) where {DimsType<:AbstractDimensions} =
-    QobjEvo(MatrixOperator(-1im * H.data.A), dims = H.dimensions, type = Operator())
-_sesolve_make_U_QobjEvo(H::QuantumObject) =
-    QobjEvo(MatrixOperator(-1im * H.data), dims = H.dimensions, type = Operator())
-_sesolve_make_U_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}) = QobjEvo(H, -1im)
+_sesolve_make_U_QobjEvo(H) = -1im * QuantumObjectEvolution(H, type = Operator())
 
 @doc raw"""
     sesolveProblem(
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 28625e6c4..5210646b2 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -113,7 +113,7 @@ function ssesolveProblem(
         sc_ops_evo_data,
     )
 
-    K = get_data(QobjEvo(H_eff_evo, -1im)) + K_l
+    K = get_data(-1im * QobjEvo(H_eff_evo)) + K_l
 
     D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
     D = DiffusionOperator(D_l)

From 34c676d849d731b42acf334344a24719afd631e2 Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Sat, 27 Sep 2025 17:20:23 +0200
Subject: [PATCH 2/5] Always update env before runtests

---
 Makefile         | 4 ++--
 test/runtests.jl | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/Makefile b/Makefile
index a7d4b1c6d..60a1f8050 100644
--- a/Makefile
+++ b/Makefile
@@ -12,10 +12,10 @@ changelog:
 	${JULIA} -e 'using Changelog; Changelog.generate(Changelog.CommonMark(), "CHANGELOG.md"; repo = "qutip/QuantumToolbox.jl")'
 
 test:
-	${JULIA} --project -e 'using Pkg; Pkg.resolve(); Pkg.test()'
+	${JULIA} --project -e 'using Pkg; Pkg.update(); Pkg.test()'
 
 docs:
-	${JULIA} --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+	${JULIA} --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.update()'
 	${JULIA} --project=docs docs/make.jl
 
 vitepress:
diff --git a/test/runtests.jl b/test/runtests.jl
index eab6514c8..96c52f274 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -26,7 +26,7 @@ const testdir = dirname(@__FILE__)
 if (GROUP == "All") || (GROUP == "Code-Quality")
     Pkg.activate("core-test/code-quality")
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
-    Pkg.instantiate()
+    Pkg.update()
 
     using QuantumToolbox
     using Aqua, JET
@@ -39,7 +39,7 @@ end
 if (GROUP == "AutoDiff_Ext")
     Pkg.activate("ext-test/cpu/autodiff")
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
-    Pkg.instantiate()
+    Pkg.update()
 
     using QuantumToolbox
     using ForwardDiff
@@ -55,7 +55,7 @@ end
 if (GROUP == "Makie_Ext")
     Pkg.activate("ext-test/cpu/makie")
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
-    Pkg.instantiate()
+    Pkg.update()
 
     using QuantumToolbox
     QuantumToolbox.about()
@@ -67,7 +67,7 @@ end
 if (GROUP == "CUDA_Ext")
     Pkg.activate("ext-test/gpu")
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
-    Pkg.instantiate()
+    Pkg.update()
 
     using QuantumToolbox
     import LinearAlgebra: Diagonal

From f7695e9bdab64978ac34870da383ed4e849f16ae Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Sat, 27 Sep 2025 18:35:24 +0200
Subject: [PATCH 3/5] First attempt to change QobjEvo generation function

---
 src/qobj/quantum_object_evo.jl | 176 ++++++++++++++++++++++-----------
 1 file changed, 120 insertions(+), 56 deletions(-)

diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index a7fd3f551..f018f8a38 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -349,72 +349,135 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
 # Arguments
 - `op_func_list::Tuple`: A tuple of tuples or operators.
 =#
-@generated function _QobjEvo_generate_data(op_func_list::Tuple)
-    op_func_list_types = op_func_list.parameters
-    N = length(op_func_list_types)
-
-    dims_expr = ()
-    func_methods_expr = ()
-    first_op = nothing
-    data_expr = :(0)
-    qobj_expr_const = :(0)
-
-    for i in 1:N
-        op_func_type = op_func_list_types[i]
-        if op_func_type <: Tuple
-            # check the structure of the tuple
-            length(op_func_type.parameters) == 2 || throw(ArgumentError("The tuple must have two elements."))
-            op_type = op_func_type.parameters[1]
-            func_type = op_func_type.parameters[2]
-            ((isoper(op_type) || issuper(op_type)) && func_type <: Function) || throw(
-                ArgumentError(
-                    "The first element must be a Operator or SuperOperator, and the second element must be a function.",
-                ),
-            )
-
-            op = :(op_func_list[$i][1])
-            dims_expr = (dims_expr..., :($op.dimensions))
-            func_methods_expr = (func_methods_expr..., :(methods(op_func_list[$i][2], [Any, Real]).ms)) # [Any, Real] means each func must accept 2 arguments
-            if i == 1
-                first_op = :($op)
-            end
-            data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i]))
+# @generated function _QobjEvo_generate_data(op_func_list::Tuple)
+#     op_func_list_types = op_func_list.parameters
+#     N = length(op_func_list_types)
+
+#     dims_expr = ()
+#     func_methods_expr = ()
+#     first_op = nothing
+#     data_expr = :(0)
+#     qobj_expr_const = :(0)
+
+#     for i in 1:N
+#         op_func_type = op_func_list_types[i]
+#         if op_func_type <: Tuple
+#             # check the structure of the tuple
+#             length(op_func_type.parameters) == 2 || throw(ArgumentError("The tuple must have two elements."))
+#             op_type = op_func_type.parameters[1]
+#             func_type = op_func_type.parameters[2]
+#             ((isoper(op_type) || issuper(op_type)) && func_type <: Function) || throw(
+#                 ArgumentError(
+#                     "The first element must be a Operator or SuperOperator, and the second element must be a function.",
+#                 ),
+#             )
+
+#             op = :(op_func_list[$i][1])
+#             dims_expr = (dims_expr..., :($op.dimensions))
+#             func_methods_expr = (func_methods_expr..., :(methods(op_func_list[$i][2], [Any, Real]).ms)) # [Any, Real] means each func must accept 2 arguments
+#             if i == 1
+#                 first_op = :($op)
+#             end
+#             data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i]))
+#         else
+#             op_type = op_func_type
+#             (isoper(op_type) || issuper(op_type)) ||
+#                 throw(ArgumentError("The element must be a Operator or SuperOperator."))
+
+#             dims_expr = (dims_expr..., :(op_func_list[$i].dimensions))
+#             if i == 1
+#                 first_op = :(op_func_list[$i])
+#             end
+#             qobj_expr_const = :($qobj_expr_const + op_func_list[$i])
+#         end
+#     end
+
+#     quote
+#         # check the dims of the operators
+#         dims = tuple($(dims_expr...))
+#         allequal(dims) || throw(ArgumentError("The dimensions of the operators must be the same."))
+
+#         # check if each func accepts 2 arguments
+#         func_methods = tuple($(func_methods_expr...))
+#         for i in eachindex(func_methods)
+#             length(func_methods[i]) == 0 && throw(
+#                 ArgumentError(
+#                     "The following function must only accept two arguments: `$(nameof(op_func_list[i][2]))(p, t)` with t<:Real",
+#                 ),
+#             )
+#         end
+
+#         data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const)
+
+#         data_expr = data_expr_const + $data_expr
+
+#         return $first_op, data_expr
+#     end
+# end
+
+Base.@constprop :aggressive function _QobjEvo_generate_data(op_func_list::Tuple)
+    first_op = _QobjEvo_get_first_op(op_func_list[1])
+
+    ops_constant = filter(op_func_list) do x
+        if x isa AbstractQuantumObject
+            x.type == first_op.type || throw(ArgumentError("The types of the operators must be the same."))
+            x.dimensions == first_op.dimensions ||
+                throw(ArgumentError("The dimensions of the operators must be the same."))
+            return true
+        elseif x isa Tuple
+            return false
         else
-            op_type = op_func_type
-            (isoper(op_type) || issuper(op_type)) ||
-                throw(ArgumentError("The element must be a Operator or SuperOperator."))
-
-            dims_expr = (dims_expr..., :(op_func_list[$i].dimensions))
-            if i == 1
-                first_op = :(op_func_list[$i])
-            end
-            qobj_expr_const = :($qobj_expr_const + op_func_list[$i])
+            throw(ArgumentError("Each element of the tuple must be either a QuantumObject or a Tuple."))
         end
     end
 
-    quote
-        # check the dims of the operators
-        dims = tuple($(dims_expr...))
-        allequal(dims) || throw(ArgumentError("The dimensions of the operators must be the same."))
-
-        # check if each func accepts 2 arguments
-        func_methods = tuple($(func_methods_expr...))
-        for i in eachindex(func_methods)
-            length(func_methods[i]) == 0 && throw(
-                ArgumentError(
-                    "The following function must only accept two arguments: `$(nameof(op_func_list[i][2]))(p, t)` with t<:Real",
-                ),
-            )
+    ops_time_dep = filter(op_func_list) do x
+        if x isa Tuple
+            _QobjEvo_check_op_func(x)
+            op = x[1]
+            _QobjEvo_check_op(op)
+            op.type == first_op.type || throw(ArgumentError("The types of the operators must be the same."))
+            op.dimensions == first_op.dimensions ||
+                throw(ArgumentError("The dimensions of the operators must be the same."))
+            return true
+        elseif x isa AbstractQuantumObject
+            return false
+        else
+            throw(ArgumentError("Each element of the tuple must be either a QuantumObject or a Tuple."))
         end
+    end
 
-        data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const)
+    data = _make_SciMLOperator(sum(ops_constant; init = zero(eltype(first_op)))) + sum(_make_SciMLOperator, ops_time_dep; init = zero(eltype(first_op)))
 
-        data_expr = data_expr_const + $data_expr
+    return first_op, data
+end
 
-        return $first_op, data_expr
-    end
+function _QobjEvo_check_op(op::AbstractQuantumObject)
+    (isoper(op) || issuper(op)) || throw(ArgumentError("The element must be a Operator or SuperOperator."))
+    return nothing
 end
 
+function _QobjEvo_check_op_func(op_func::Tuple)
+    length(op_func) == 2 || throw(ArgumentError("The tuple must have two elements."))
+    _QobjEvo_check_op(op_func[1])
+    (op_func[2] isa Function) || throw(ArgumentError("The second element must be a function."))
+    methods(op_func[2], (Any, Real)) |> length == 0 &&
+        throw(ArgumentError("The following function must only accept two arguments: `$(nameof(op_func[2]))(p, t)` with t<:Real"))
+    return nothing
+end
+
+_QobjEvo_get_first_op(op_func_list_1::Union{Tuple,AbstractQuantumObject}) =
+    if op_func_list_1 isa Tuple
+        _QobjEvo_check_op_func(op_func_list_1)
+        op = op_func_list_1[1]
+        _QobjEvo_check_op(op)
+        return op
+    else
+        op = op_func_list_1
+        _QobjEvo_check_op(op)
+        return op
+    end
+
 function _make_SciMLOperator(op_func::Tuple)
     T = eltype(op_func[1])
     update_func = (a, u, p, t) -> op_func[2](p, t)
@@ -422,6 +485,7 @@ function _make_SciMLOperator(op_func::Tuple)
 end
 
 _make_SciMLOperator(op::AbstractQuantumObject) = _promote_to_scimloperator(op.data)
+_make_SciMLOperator(op::Number) = op
 
 _promote_to_scimloperator(data::AbstractMatrix) = MatrixOperator(data)
 _promote_to_scimloperator(data::AbstractSciMLOperator) = data

From 371af73423cc7bbe292419f48f631938eab43d7d Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Sat, 27 Sep 2025 20:01:15 +0200
Subject: [PATCH 4/5] Make QobjEvo generation working

---
 src/qobj/quantum_object_evo.jl | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index f018f8a38..5725445bf 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -447,7 +447,11 @@ Base.@constprop :aggressive function _QobjEvo_generate_data(op_func_list::Tuple)
         end
     end
 
-    data = _make_SciMLOperator(sum(ops_constant; init = zero(eltype(first_op)))) + sum(_make_SciMLOperator, ops_time_dep; init = zero(eltype(first_op)))
+    data_const = sum(_make_SciMLOperator, ops_constant; init = zero(eltype(first_op))*I)
+    data_td =
+        length(ops_time_dep) == 1 ? _make_SciMLOperator(ops_time_dep[1]) :
+        AddedOperator(map(_make_SciMLOperator, ops_time_dep))
+    data = isempty(ops_constant) ? data_td : data_const + data_td
 
     return first_op, data
 end
@@ -461,8 +465,11 @@ function _QobjEvo_check_op_func(op_func::Tuple)
     length(op_func) == 2 || throw(ArgumentError("The tuple must have two elements."))
     _QobjEvo_check_op(op_func[1])
     (op_func[2] isa Function) || throw(ArgumentError("The second element must be a function."))
-    methods(op_func[2], (Any, Real)) |> length == 0 &&
-        throw(ArgumentError("The following function must only accept two arguments: `$(nameof(op_func[2]))(p, t)` with t<:Real"))
+    methods(op_func[2], (Any, Real)) |> length == 0 && throw(
+        ArgumentError(
+            "The following function must only accept two arguments: `$(nameof(op_func[2]))(p, t)` with t<:Real",
+        ),
+    )
     return nothing
 end
 
@@ -479,13 +486,12 @@ _QobjEvo_get_first_op(op_func_list_1::Union{Tuple,AbstractQuantumObject}) =
     end
 
 function _make_SciMLOperator(op_func::Tuple)
-    T = eltype(op_func[1])
-    update_func = (a, u, p, t) -> op_func[2](p, t)
-    return ScalarOperator(zero(T), update_func) * _promote_to_scimloperator(op_func[1].data)
+    op, coef = op_func
+    T = eltype(op)
+    update_func = (a, u, p, t) -> coef(p, t)
+    return ScalarOperator(zero(T), update_func) * _promote_to_scimloperator(op.data)
 end
-
 _make_SciMLOperator(op::AbstractQuantumObject) = _promote_to_scimloperator(op.data)
-_make_SciMLOperator(op::Number) = op
 
 _promote_to_scimloperator(data::AbstractMatrix) = MatrixOperator(data)
 _promote_to_scimloperator(data::AbstractSciMLOperator) = data

From 0b7a2b5c65c5eb9f579924425e19432ed9cd8abc Mon Sep 17 00:00:00 2001
From: Alberto Mercurio <alberto.mercurio96@gmail.com>
Date: Sat, 27 Sep 2025 21:21:17 +0200
Subject: [PATCH 5/5] Fix type-instabilities error

---
 src/qobj/quantum_object_evo.jl | 76 +++-------------------------------
 src/time_evolution/mcsolve.jl  | 10 ++---
 src/time_evolution/smesolve.jl |  2 +-
 src/time_evolution/ssesolve.jl |  2 +-
 4 files changed, 12 insertions(+), 78 deletions(-)

diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index 5725445bf..8e4de483b 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -349,77 +349,11 @@ Parse the `op_func_list` and generate the data for the `QuantumObjectEvolution`
 # Arguments
 - `op_func_list::Tuple`: A tuple of tuples or operators.
 =#
-# @generated function _QobjEvo_generate_data(op_func_list::Tuple)
-#     op_func_list_types = op_func_list.parameters
-#     N = length(op_func_list_types)
-
-#     dims_expr = ()
-#     func_methods_expr = ()
-#     first_op = nothing
-#     data_expr = :(0)
-#     qobj_expr_const = :(0)
-
-#     for i in 1:N
-#         op_func_type = op_func_list_types[i]
-#         if op_func_type <: Tuple
-#             # check the structure of the tuple
-#             length(op_func_type.parameters) == 2 || throw(ArgumentError("The tuple must have two elements."))
-#             op_type = op_func_type.parameters[1]
-#             func_type = op_func_type.parameters[2]
-#             ((isoper(op_type) || issuper(op_type)) && func_type <: Function) || throw(
-#                 ArgumentError(
-#                     "The first element must be a Operator or SuperOperator, and the second element must be a function.",
-#                 ),
-#             )
-
-#             op = :(op_func_list[$i][1])
-#             dims_expr = (dims_expr..., :($op.dimensions))
-#             func_methods_expr = (func_methods_expr..., :(methods(op_func_list[$i][2], [Any, Real]).ms)) # [Any, Real] means each func must accept 2 arguments
-#             if i == 1
-#                 first_op = :($op)
-#             end
-#             data_expr = :($data_expr + _make_SciMLOperator(op_func_list[$i]))
-#         else
-#             op_type = op_func_type
-#             (isoper(op_type) || issuper(op_type)) ||
-#                 throw(ArgumentError("The element must be a Operator or SuperOperator."))
-
-#             dims_expr = (dims_expr..., :(op_func_list[$i].dimensions))
-#             if i == 1
-#                 first_op = :(op_func_list[$i])
-#             end
-#             qobj_expr_const = :($qobj_expr_const + op_func_list[$i])
-#         end
-#     end
-
-#     quote
-#         # check the dims of the operators
-#         dims = tuple($(dims_expr...))
-#         allequal(dims) || throw(ArgumentError("The dimensions of the operators must be the same."))
-
-#         # check if each func accepts 2 arguments
-#         func_methods = tuple($(func_methods_expr...))
-#         for i in eachindex(func_methods)
-#             length(func_methods[i]) == 0 && throw(
-#                 ArgumentError(
-#                     "The following function must only accept two arguments: `$(nameof(op_func_list[i][2]))(p, t)` with t<:Real",
-#                 ),
-#             )
-#         end
-
-#         data_expr_const = $qobj_expr_const isa Integer ? $qobj_expr_const : _make_SciMLOperator($qobj_expr_const)
-
-#         data_expr = data_expr_const + $data_expr
-
-#         return $first_op, data_expr
-#     end
-# end
-
-Base.@constprop :aggressive function _QobjEvo_generate_data(op_func_list::Tuple)
+function _QobjEvo_generate_data(op_func_list::Tuple)
     first_op = _QobjEvo_get_first_op(op_func_list[1])
 
     ops_constant = filter(op_func_list) do x
-        if x isa AbstractQuantumObject
+        if x isa QuantumObject
             x.type == first_op.type || throw(ArgumentError("The types of the operators must be the same."))
             x.dimensions == first_op.dimensions ||
                 throw(ArgumentError("The dimensions of the operators must be the same."))
@@ -440,18 +374,18 @@ Base.@constprop :aggressive function _QobjEvo_generate_data(op_func_list::Tuple)
             op.dimensions == first_op.dimensions ||
                 throw(ArgumentError("The dimensions of the operators must be the same."))
             return true
-        elseif x isa AbstractQuantumObject
+        elseif x isa QuantumObject
             return false
         else
             throw(ArgumentError("Each element of the tuple must be either a QuantumObject or a Tuple."))
         end
     end
 
-    data_const = sum(_make_SciMLOperator, ops_constant; init = zero(eltype(first_op))*I)
+    data_const = isempty(ops_constant) ? zero(eltype(first_op)) : _make_SciMLOperator(sum(ops_constant))
     data_td =
         length(ops_time_dep) == 1 ? _make_SciMLOperator(ops_time_dep[1]) :
         AddedOperator(map(_make_SciMLOperator, ops_time_dep))
-    data = isempty(ops_constant) ? data_td : data_const + data_td
+    data = data_const + data_td
 
     return first_op, data
 end
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 98b60076e..ef195b591 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -26,16 +26,16 @@ function _normalize_state!(u, dims, normalize_states)
 end
 
 function _mcsolve_make_Heff_QobjEvo(H::QuantumObject, c_ops)
-    c_ops isa Nothing && return QobjEvo(H)
-    return QobjEvo(H - 1im * mapreduce(op -> op' * op, +, c_ops) / 2)
+    c_ops isa Nothing && return QuantumObjectEvolution(H)
+    return QuantumObjectEvolution(H - 1im * mapreduce(op -> op' * op, +, c_ops) / 2)
 end
 function _mcsolve_make_Heff_QobjEvo(H::Tuple, c_ops)
-    c_ops isa Nothing && return QobjEvo(H)
-    return QobjEvo((H..., -1im * sum(op -> op' * op, c_ops) / 2))
+    c_ops isa Nothing && return QuantumObjectEvolution(H)
+    return QuantumObjectEvolution((H..., -1im * sum(op -> op' * op, c_ops) / 2))
 end
 function _mcsolve_make_Heff_QobjEvo(H::QuantumObjectEvolution, c_ops)
     c_ops isa Nothing && return H
-    return H + QobjEvo(sum(op -> -1im * op' * op / 2, c_ops))
+    return H + QuantumObjectEvolution(sum(op -> -1im * op' * op / 2, c_ops))
 end
 
 @doc raw"""
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index a09ba9437..2c969f7bc 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -119,7 +119,7 @@ function smesolveProblem(
         # TODO: # Currently, we are assuming a time-independent MatrixOperator
         # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2
         op_vec = mat2vec(adjoint(op.A))
-        return _spre(op, Id) + _spost(op', Id) + _smesolve_ScalarOperator(op_vec) * Id_op
+        return AddedOperator(_spre(op, Id), _spost(op', Id), _smesolve_ScalarOperator(op_vec) * Id_op)
     end
     D = DiffusionOperator(D_l)
 
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 5210646b2..99c143795 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -113,7 +113,7 @@ function ssesolveProblem(
         sc_ops_evo_data,
     )
 
-    K = get_data(-1im * QobjEvo(H_eff_evo)) + K_l
+    K = get_data(-1im * QuantumObjectEvolution(H_eff_evo)) + K_l
 
     D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
     D = DiffusionOperator(D_l)