Merge branch 'main' into fix/makie-warning

ytdHuang · ytdHuang · commit 5dc1a9b34ab5 · 2025-09-27T17:59:31.000+08:00
diff --git a/.github/workflows/Benchmarks.yml b/.github/workflows/Benchmarks.yml
@@ -34,7 +34,7 @@ jobs:
     runs-on: ubuntu-latest
     if: ${{ !github.event.pull_request.draft }}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
diff --git a/.github/workflows/CI-Julia-nightly.yml b/.github/workflows/CI-Julia-nightly.yml
@@ -48,7 +48,7 @@ jobs:
           - 'Core'
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -75,7 +75,7 @@ jobs:
             group: 'AutoDiff_Ext'
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
diff --git a/.github/workflows/ChangeLogCheck.yml b/.github/workflows/ChangeLogCheck.yml
@@ -16,7 +16,7 @@ jobs:
           skipLabels: 'Skip ChangeLog'
       
       # check whether the format of CHANGELOG.md is correct
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
diff --git a/.github/workflows/CleanPreviewDoc.yml b/.github/workflows/CleanPreviewDoc.yml
@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout gh-pages branch
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
         with:
           ref: gh-pages
       - name: Delete preview and history + push changes
diff --git a/.github/workflows/Code-Quality.yml b/.github/workflows/Code-Quality.yml
@@ -48,7 +48,7 @@ jobs:
           - 'Code-Quality'
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
@@ -25,7 +25,7 @@ jobs:
     runs-on: ubuntu-latest
     if: ${{ !github.event.pull_request.draft }}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
@@ -8,6 +8,6 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout Actions Repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v5
       - name: Check spelling
-        uses: crate-ci/typos@v1.35.3
+        uses: crate-ci/typos@v1.36.2
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
@@ -43,7 +43,7 @@ jobs:
     runs-on: ubuntu-latest
     if: ${{ !github.event.pull_request.draft }}
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,7 +7,20 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+## [v0.35.0]
+Release date: 2025-09-03
+
+- Add support of `QobjEvo` for `steadystate` (ODE solver only). ([#536])
+- Changes to `SteadyStateODESolver`. ([#537])
+  - Introduce the tolerances for `steadystate` terminate condition (two new fields: `terminate_reltol = 1e-5` and `terminate_abstol = 1e-7`)
+  - Fix keyword argument handling for `SteadyStateODESolver` before passing to `mesolve`.
+- Fix incorrect `negativity` and `partial_transpose` for arbitrary subsystem dimension. ([#539])
+
+## [v0.34.1]
+Release date: 2025-08-23
+
 - Improve Bloch sphere rendering for animation. ([#520])
+- Add support to `Enzyme.jl` for `sesolve` and `mesolve`. ([#531])
 
 ## [v0.34.0]
 Release date: 2025-07-29
@@ -208,6 +221,8 @@ Release date: 2024-11-13
 [v0.32.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.32.1
 [v0.33.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.33.0
 [v0.34.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.34.0
+[v0.34.1]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.34.1
+[v0.35.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.35.0
 [#86]: https://github.com/qutip/QuantumToolbox.jl/issues/86
 [#139]: https://github.com/qutip/QuantumToolbox.jl/issues/139
 [#271]: https://github.com/qutip/QuantumToolbox.jl/issues/271
@@ -295,3 +310,7 @@ Release date: 2024-11-13
 [#515]: https://github.com/qutip/QuantumToolbox.jl/issues/515
 [#517]: https://github.com/qutip/QuantumToolbox.jl/issues/517
 [#520]: https://github.com/qutip/QuantumToolbox.jl/issues/520
+[#531]: https://github.com/qutip/QuantumToolbox.jl/issues/531
+[#536]: https://github.com/qutip/QuantumToolbox.jl/issues/536
+[#537]: https://github.com/qutip/QuantumToolbox.jl/issues/537
+[#539]: https://github.com/qutip/QuantumToolbox.jl/issues/539
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumToolbox"
 uuid = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 authors = ["Alberto Mercurio", "Yi-Te Huang"]
-version = "0.34.0"
+version = "0.35.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
diff --git a/src/negativity.jl b/src/negativity.jl
@@ -85,20 +85,20 @@ end
 
 # for dense matrices
 function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
-    mask2 = [1 + Int(i) for i in mask]
+    nsys = length(mask)
+    mask2 = reverse([mask[s] ? 2 : 1 for s in 1:nsys])
     # mask2 has elements with values equal to 1 or 2
-    #   1 - the subsystem don't need to be transposed
-    #   2 - the subsystem need be transposed
+    #   1 - the subsystem (in reversed order) don't need to be transposed
+    #   2 - the subsystem (in reversed order) need to be transposed
 
-    nsys = length(mask2)
-    dims = dimensions_to_dims(get_dimensions_to(ρ))
+    dims_rev = reverse(dimensions_to_dims(get_dimensions_to(ρ)))
     pt_dims = reshape(Vector(1:(2*nsys)), (nsys, 2))
     pt_idx = [
         [pt_dims[n, mask2[n]] for n in 1:nsys]   # origin   value in mask2
         [pt_dims[n, 3-mask2[n]] for n in 1:nsys] # opposite value in mask2 (1 -> 2, and 2 -> 1)
     ]
     return QuantumObject(
-        reshape(permutedims(reshape(ρ.data, (dims..., dims...)), pt_idx), size(ρ)),
+        reshape(permutedims(reshape(ρ.data, (dims_rev..., dims_rev...)), pt_idx), size(ρ)),
         Operator(),
         Dimensions(ρ.dimensions.to),
     )
@@ -110,7 +110,8 @@ function _partial_transpose(
     mask::Vector{Bool},
 ) where {DimsType<:AbstractDimensions}
     M, N = size(ρ)
-    dimsTuple = Tuple(dimensions_to_dims(get_dimensions_to(ρ)))
+    dims_rev = reverse(Tuple(dimensions_to_dims(get_dimensions_to(ρ))))
+    mask_rev = reverse(mask)
     colptr = ρ.data.colptr
     rowval = ρ.data.rowval
     nzval = ρ.data.nzval
@@ -130,13 +131,13 @@ function _partial_transpose(
                 I_pt[n] = i
                 J_pt[n] = j
             else
-                ket_pt = [Base._ind2sub(dimsTuple, i)...]
-                bra_pt = [Base._ind2sub(dimsTuple, j)...]
-                for sys in findall(m -> m, mask)
+                ket_pt = [Base._ind2sub(dims_rev, i)...]
+                bra_pt = [Base._ind2sub(dims_rev, j)...]
+                for sys in findall(m -> m, mask_rev)
                     @inbounds ket_pt[sys], bra_pt[sys] = bra_pt[sys], ket_pt[sys]
                 end
-                I_pt[n] = Base._sub2ind(dimsTuple, ket_pt...)
-                J_pt[n] = Base._sub2ind(dimsTuple, bra_pt...)
+                I_pt[n] = Base._sub2ind(dims_rev, ket_pt...)
+                J_pt[n] = Base._sub2ind(dims_rev, bra_pt...)
             end
             V_pt[n] = nzval[p]
         end
diff --git a/src/steadystate.jl b/src/steadystate.jl
@@ -47,11 +47,11 @@ end
         alg = Tsit5(),
         ψ0 = nothing,
         tmax = Inf,
-        )
-
-An ordinary differential equation (ODE) solver for solving [`steadystate`](@ref).
+        terminate_reltol = 1e-5,
+        terminate_abstol = 1e-7
+    )
 
-Solve the stationary state based on time evolution (ordinary differential equations; `OrdinaryDiffEq.jl`) with a given initial state.
+An ordinary differential equation (ODE) solver for solving [`steadystate`](@ref). It solves the stationary state based on [`mesolve`](@ref) with a termination condition.
 
 The termination condition of the stationary state ``|\rho\rangle\rangle`` is that either the following condition is `true`:
 
@@ -69,14 +69,26 @@ or
 - `alg::OrdinaryDiffEqAlgorithm=Tsit5()`: The algorithm to solve the ODE.
 - `ψ0::Union{Nothing,QuantumObject}=nothing`: The initial state of the system. If not specified, a random pure state will be generated.
 - `tmax::Real=Inf`: The final time step for the steady state problem.
+- `terminate_reltol` = The relative tolerance for stationary state terminate condition. Default to `1e-5`.
+- `terminate_abstol` = The absolute tolerance for stationary state terminate condition. Default to `1e-7`.
+
+!!! warning "Tolerances for terminate condition"
+    The terminate condition tolerances `terminate_reltol` and `terminate_abstol` should be larger than `reltol` and `abstol` of [`mesolve`](@ref), respectively.
 
-For more details about the solvers, please refer to [`OrdinaryDiffEq.jl`](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
+For more details about the solving `alg`orithms, please refer to [`OrdinaryDiffEq.jl`](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 """
-Base.@kwdef struct SteadyStateODESolver{MT<:OrdinaryDiffEqAlgorithm,ST<:Union{Nothing,QuantumObject},T<:Real} <:
-                   SteadyStateSolver
+Base.@kwdef struct SteadyStateODESolver{
+    MT<:OrdinaryDiffEqAlgorithm,
+    ST<:Union{Nothing,QuantumObject},
+    TT<:Real,
+    RT<:Real,
+    AT<:Real,
+} <: SteadyStateSolver
     alg::MT = Tsit5()
     ψ0::ST = nothing
-    tmax::T = Inf
+    tmax::TT = Inf
+    terminate_reltol::RT = 10 * DEFAULT_ODE_SOLVER_OPTIONS.reltol
+    terminate_abstol::AT = 10 * DEFAULT_ODE_SOLVER_OPTIONS.abstol
 end
 
 @doc raw"""
@@ -96,7 +108,7 @@ end
 
 @doc raw"""
     steadystate(
-        H::QuantumObject{OpType},
+        H::AbstractQuantumObject{OpType},
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateSolver = SteadyStateDirectSolver(),
         kwargs...,
@@ -111,7 +123,7 @@ Solve the stationary state based on different solvers.
 - `kwargs`: The keyword arguments for the solver.
 """
 function steadystate(
-    H::QuantumObject{OpType},
+    H::AbstractQuantumObject{OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
@@ -199,30 +211,40 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateDirect
     return QuantumObject(ρss, Operator(), L.dimensions)
 end
 
-function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...)
-    tmax = solver.tmax
-
+function _steadystate(L::AbstractQuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...)
     ψ0 = isnothing(solver.ψ0) ? rand_ket(L.dimensions) : solver.ψ0
-    abstol = haskey(kwargs, :abstol) ? kwargs[:abstol] : DEFAULT_ODE_SOLVER_OPTIONS.abstol
-    reltol = haskey(kwargs, :reltol) ? kwargs[:reltol] : DEFAULT_ODE_SOLVER_OPTIONS.reltol
-
     ftype = _float_type(ψ0)
-    _terminate_func = SteadyStateODECondition(similar(mat2vec(ket2dm(ψ0)).data))
-    cb = TerminateSteadyState(abstol, reltol, _terminate_func)
-    sol = mesolve(
-        L,
-        ψ0,
-        [ftype(0), ftype(tmax)],
-        progress_bar = Val(false),
-        save_everystep = false,
-        saveat = ftype[],
-        callback = cb,
+    tlist = [ftype(0), ftype(solver.tmax)]
+
+    # overwrite some kwargs and throw warning message to tell the users that we are ignoring these settings
+    haskey(kwargs, :progress_bar) && @warn "Ignore keyword argument 'progress_bar' for SteadyStateODESolver"
+    haskey(kwargs, :save_everystep) && @warn "Ignore keyword argument 'save_everystep' for SteadyStateODESolver"
+    haskey(kwargs, :saveat) && @warn "Ignore keyword argument 'saveat' for SteadyStateODESolver"
+    kwargs2 = merge(
+        NamedTuple(kwargs), # we convert to NamedTuple just in case if kwargs is empty
+        (progress_bar = Val(false), save_everystep = false, saveat = ftype[]),
+    )
+
+    # add terminate condition (callback)
+    cb = TerminateSteadyState(
+        solver.terminate_abstol,
+        solver.terminate_reltol,
+        SteadyStateODECondition(similar(mat2vec(ket2dm(ψ0)).data)),
     )
+    kwargs3 = _merge_kwargs_with_callback(kwargs2, cb)
 
+    sol = mesolve(L, ψ0, tlist; kwargs3...)
     ρss = sol.states[end]
     return ρss
 end
 
+_steadystate(
+    L::QuantumObjectEvolution{SuperOperator},
+    solver::T;
+    kwargs...,
+) where {T<:Union{SteadyStateDirectSolver,SteadyStateEigenSolver,SteadyStateLinearSolver}} =
+    throw(ArgumentError("$(get_typename_wrapper(solver)) does not support QobjEvo."))
+
 struct SteadyStateODECondition{CT<:AbstractArray}
     cache::CT
 end
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
@@ -414,6 +414,8 @@ function mcsolve(
     col_times = map(i -> _mc_get_jump_callback(sol[:, i]).affect!.col_times, eachindex(sol))
     col_which = map(i -> _mc_get_jump_callback(sol[:, i]).affect!.col_which, eachindex(sol))
 
+    kwargs = NamedTuple(_sol_1.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility
+
     return TimeEvolutionMCSol(
         ntraj,
         ens_prob_mc.times,
@@ -424,7 +426,7 @@ function mcsolve(
         col_which,
         sol.converged,
         _sol_1.alg,
-        _sol_1.prob.kwargs[:abstol],
-        _sol_1.prob.kwargs[:reltol],
+        kwargs.abstol,
+        kwargs.reltol,
     )
 end
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
@@ -178,6 +178,11 @@ function mesolve(
         kwargs...,
     )
 
+    # Move sensealg argument to solve for Enzyme.jl support.
+    # TODO: Remove it when https://github.com/SciML/SciMLSensitivity.jl/issues/1225 is fixed.
+    sensealg = get(kwargs, :sensealg, nothing)
+    kwargs_filtered = isnothing(sensealg) ? kwargs : Base.structdiff((; kwargs...), (sensealg = sensealg,))
+
     prob = mesolveProblem(
         H,
         ψ0,
@@ -188,14 +193,19 @@ function mesolve(
         params = params,
         progress_bar = progress_bar,
         inplace = inplace,
-        kwargs...,
+        kwargs_filtered...,
     )
 
-    return mesolve(prob, alg)
+    # TODO: Remove sensealg when https://github.com/SciML/SciMLSensitivity.jl/issues/1225 is fixed
+    if isnothing(sensealg)
+        return mesolve(prob, alg)
+    else
+        return mesolve(prob, alg; sensealg = sensealg)
+    end
 end
 
-function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5())
-    sol = solve(prob.prob, alg)
+function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); kwargs...)
+    sol = solve(prob.prob, alg; kwargs...)
 
     # No type instabilities since `isoperket` is a Val, and so it is known at compile time
     if getVal(prob.kwargs.isoperket)
@@ -204,14 +214,16 @@ function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
         ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator(), dims = prob.dimensions), sol.u)
     end
 
+    kwargs = NamedTuple(sol.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility
+
     return TimeEvolutionSol(
         prob.times,
         sol.t,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
         sol.alg,
-        sol.prob.kwargs[:abstol],
-        sol.prob.kwargs[:reltol],
+        kwargs.abstol,
+        kwargs.reltol,
     )
 end
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
diff --git a/test/core-test/negativity_and_partial_transpose.jl b/test/core-test/negativity_and_partial_transpose.jl
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl
