Merge branch 'main' into julia-1.12

Author: ytdHuang

.github/workflows/Benchmarks.yml

@@ -37,7 +37,7 @@ jobs:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1'
+          version: '1.11'
           arch: x64
       - uses: julia-actions/cache@v2
       - name: Run benchmark

.github/workflows/CleanPreviewDoc.yml

@@ -15,7 +15,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(s) and reset history

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v5
       - name: Check spelling
-        uses: crate-ci/typos@v1.37.2
+        uses: crate-ci/typos@v1.38.1

CHANGELOG.md

@@ -7,7 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Introduce new methods of `sesolve_map` and `mesolve_map` for advanced usage. Users can now customize their own `iter`ator structure, `prob_func` and `output_func`. ([#565])
+
+## [v0.37.0]
+Release date: 2025-10-12
+
 - Fix `cite()` bibtex output. ([#552])
+- Implement `sesolve_map` and `mesolve_map` for solving multiple initial states and parameter sets in parallel. ([#554])
 - Add `qeye_like` and `qzero_like`, which are synonyms of `one` and `zero`. ([#555])
 - Add steadystate and DSF benchmarks. The `SteadyStateODESOlver` tolerances are lowered to `terminate_reltol=1e-4` and `terminate_abstol=1e-6` to improve speed at the cost of accuracy. ([#557])
 
@@ -234,6 +240,7 @@ Release date: 2024-11-13
 [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
 [v0.36.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.36.0
+[v0.37.0]: https://github.com/qutip/QuantumToolbox.jl/releases/tag/v0.37.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
@@ -328,5 +335,7 @@ Release date: 2024-11-13
 [#544]: https://github.com/qutip/QuantumToolbox.jl/issues/544
 [#546]: https://github.com/qutip/QuantumToolbox.jl/issues/546
 [#552]: https://github.com/qutip/QuantumToolbox.jl/issues/552
+[#554]: https://github.com/qutip/QuantumToolbox.jl/issues/554
 [#555]: https://github.com/qutip/QuantumToolbox.jl/issues/555
 [#557]: https://github.com/qutip/QuantumToolbox.jl/issues/557
+[#565]: https://github.com/qutip/QuantumToolbox.jl/issues/565

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumToolbox"
 uuid = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 authors = ["Alberto Mercurio", "Yi-Te Huang"]
-version = "0.36.0"
+version = "0.37.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

benchmarks/Project.toml

@@ -1,4 +1,8 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

benchmarks/autodiff.jl

@@ -0,0 +1,92 @@
+function benchmark_autodiff!(SUITE)
+    # Use harmonic oscillator system for both sesolve and mesolve
+    N = 20
+    a = destroy(N)
+    ψ0 = fock(N, 0)
+    tlist = range(0, 40, 100)
+
+    # ---- SESOLVE ----
+    # For direct Forward differentiation
+    function my_f_sesolve_direct(p)
+        H = p[1] * a' * a + p[2] * (a + a')
+        sol = sesolve(H, ψ0, tlist, progress_bar = Val(false))
+        return real(expect(a' * a, sol.states[end]))
+    end
+
+    # For SciMLSensitivity.jl (reverse mode with Zygote and Enzyme)
+    coef_Δ(p, t) = p[1]
+    coef_F(p, t) = p[2]
+    H_evo = QobjEvo(a' * a, coef_Δ) + QobjEvo(a + a', coef_F)
+
+    function my_f_sesolve(p)
+        sol = sesolve(
+            H_evo,
+            ψ0,
+            tlist,
+            progress_bar = Val(false),
+            params = p,
+            sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()),
+        )
+        return real(expect(a' * a, sol.states[end]))
+    end
+
+    # ---- MESOLVE ----
+    # For direct Forward differentiation
+    function my_f_mesolve_direct(p)
+        H = p[1] * a' * a + p[2] * (a + a')
+        c_ops = [sqrt(p[3]) * a]
+        sol = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
+        return real(expect(a' * a, sol.states[end]))
+    end
+
+    # For SciMLSensitivity.jl (reverse mode with Zygote and Enzyme)
+    coef_γ(p, t) = sqrt(p[3])
+    c_ops = [QobjEvo(a, coef_γ)]
+    L = liouvillian(H_evo, c_ops)
+
+    function my_f_mesolve(p)
+        sol = mesolve(
+            L,
+            ψ0,
+            tlist,
+            progress_bar = Val(false),
+            params = p,
+            sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()),
+        )
+        return real(expect(a' * a, sol.states[end]))
+    end
+
+    # Parameters for benchmarks
+    params_sesolve = [1.0, 1.0]
+    params_mesolve = [1.0, 1.0, 1.0]
+
+    # Benchmark sesolve - Forward
+    SUITE["Autodiff"]["sesolve"]["Forward"] = @benchmarkable ForwardDiff.gradient($my_f_sesolve_direct, $params_sesolve)
+
+    # Benchmark sesolve - Reverse (Zygote)
+    SUITE["Autodiff"]["sesolve"]["Reverse (Zygote)"] = @benchmarkable Zygote.gradient($my_f_sesolve, $params_sesolve)
+
+    # Benchmark sesolve - Reverse (Enzyme)
+    SUITE["Autodiff"]["sesolve"]["Reverse (Enzyme)"] = @benchmarkable Enzyme.autodiff(
+        Enzyme.set_runtime_activity(Enzyme.Reverse),
+        Const($my_f_sesolve),
+        Active,
+        Duplicated($params_sesolve, dparams_sesolve),
+    ) setup=(dparams_sesolve = Enzyme.make_zero($params_sesolve))
+
+    # Benchmark mesolve - Forward
+    SUITE["Autodiff"]["mesolve"]["Forward"] = @benchmarkable ForwardDiff.gradient($my_f_mesolve_direct, $params_mesolve)
+
+    # Benchmark mesolve - Reverse (Zygote)
+    SUITE["Autodiff"]["mesolve"]["Reverse (Zygote)"] = @benchmarkable Zygote.gradient($my_f_mesolve, $params_mesolve)
+
+    # Benchmark mesolve - Reverse (Enzyme)
+    SUITE["Autodiff"]["mesolve"]["Reverse (Enzyme)"] = @benchmarkable Enzyme.autodiff(
+        Enzyme.set_runtime_activity(Enzyme.Reverse),
+        Const($my_f_mesolve),
+        Active,
+        Duplicated($params_mesolve, dparams_mesolve),
+    ) setup=(dparams_mesolve = Enzyme.make_zero($params_mesolve))
+
+    return nothing
+end

benchmarks/runbenchmarks.jl

@@ -3,6 +3,10 @@ using LinearAlgebra
 using SparseArrays
 using QuantumToolbox
 using SciMLBase: EnsembleSerial, EnsembleThreads
+using ForwardDiff
+using Zygote
+using Enzyme: Enzyme, Const, Active, Duplicated
+using SciMLSensitivity: BacksolveAdjoint, EnzymeVJP
 
 BLAS.set_num_threads(1)
 
@@ -14,13 +18,15 @@ include("dynamical_shifted_fock.jl")
 include("eigenvalues.jl")
 include("steadystate.jl")
 include("timeevolution.jl")
+include("autodiff.jl")
 
 benchmark_correlations_and_spectrum!(SUITE)
 benchmark_dfd!(SUITE)
 benchmark_dsf!(SUITE)
 benchmark_eigenvalues!(SUITE)
 benchmark_steadystate!(SUITE)
 benchmark_timeevolution!(SUITE)
+benchmark_autodiff!(SUITE)
 
 BenchmarkTools.tune!(SUITE)
 results = BenchmarkTools.run(SUITE, verbose = true)

docs/src/resources/api.md

@@ -215,6 +215,8 @@ mesolve
 mcsolve
 ssesolve
 smesolve
+sesolve_map
+mesolve_map
 dfd_mesolve
 liouvillian
 liouvillian_generalized

docs/src/users_guide/cluster.md

@@ -193,72 +193,44 @@ flush(stdout)
 end
 
 @everywhere begin
-    const Nc = 20
-    const ωc = 1.0
-    const g = 0.05
-    const γ = 0.01
-    const F = 0.01
+    Nc = 20
+    ωc = 1.0
+    g = 0.05
+    γ = 0.01
+    F = 0.01
 
-    const a = tensor(destroy(Nc), qeye(2))
+    a = tensor(destroy(Nc), qeye(2))
 
-    const σm = tensor(qeye(Nc), sigmam())
-    const σp = tensor(qeye(Nc), sigmap())
+    σm = tensor(qeye(Nc), sigmam())
+    σp = tensor(qeye(Nc), sigmap())
+    σz = tensor(qeye(Nc), sigmaz())
 
-    H(ωq) = ωc * a' * a + ωq * tensor(num(Nc), qeye(2)) + g * (a' * σm + a * σp)
+    ωq_fun(p, t) = p[1] # ωq
+    drive_fun(p, t) = p[3] * cos(p[2] * t) # F cos(ωd t)
 
-    coef(p, t) = p.F * cos(p.ωd * t) # coefficient for the time-dependent term
+    H = ωc * a' * a + QobjEvo(σz / 2, ωq_fun) + g * (a' * σm + a * σp) + QobjEvo(a + a', drive_fun)
 
-    const c_ops = [sqrt(γ) * a, sqrt(γ) * σm]
-    const e_ops = [a' * a]
+    c_ops = [sqrt(γ) * a, sqrt(γ) * σm]
+    e_ops = [a' * a]
 end
 
-# Define the ODE problem and the EnsembleProblem generation function
+ωq_list = range(ωc - 3*g, ωc + 3*g, 100)
+ωd_list = range(ωc - 3*g, ωc + 3*g, 100)
+F_list = [F]
 
-@everywhere begin
-    ωq_list = range(ωc - 3*g, ωc + 3*g, 100)
-    ωd_list = range(ωc - 3*g, ωc + 3*g, 100)
-
-    const iter = collect(Iterators.product(ωq_list, ωd_list))
-
-    function my_prob_func(prob, i, repeat, channel)
-        ωq, ωd = iter[i]
-        H_i = H(ωq)
-        H_d_i = H_i + QobjEvo(a + a', coef) # Hamiltonian with a driving term
-
-        L = liouvillian(H_d_i, c_ops).data # Make the Liouvillian
-
-        put!(channel, true) # Update the progress bar channel
-
-        remake(prob, f=L, p=(F = F, ωd = ωd))
-    end
-end
-
-ωq, ωd = iter[1]
-H0 = H(ωq) + QobjEvo(a + a', coef)
-ψ0 = tensor(fock(Nc, 0), basis(2, 1)) # Ground State
+ψ_list = [
+    tensor(fock(Nc, 0), basis(2, 0)),
+    tensor(fock(Nc, 0), basis(2, 1))
+]
 tlist = range(0, 20 / γ, 1000)
 
-prob = mesolveProblem(H0, ψ0, tlist, c_ops, e_ops=e_ops, progress_bar=Val(false), params=(F = F, ωd = ωd))
-
-### Just to print the progress bar
-progr = ProgressBar(length(iter))
-progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
-###
-ens_prob = EnsembleProblem(prob.prob, prob_func=(prob, i, repeat) -> my_prob_func(prob, i, repeat, progr_channel))
-
-
-@sync begin
-    @async while take!(progr_channel)
-        next!(progr)
-    end
-
-    @async begin
-        sol = solve(ens_prob, Tsit5(), EnsembleSplitThreads(), trajectories = length(iter))
-        put!(progr_channel, false)
-    end
-end
+sol = mesolve_map(H, ψ_list, tlist, c_ops; e_ops = e_ops, params = (ωq_list, ωd_list, F_list), ensemblealg = EnsembleSplitThreads())
 ```
 
-We are using the [`mesolveProblem`](@ref) function to define the master equation problem. We added some code to manage the progress bar, which is updated through a `RemoteChannel`. The `prob_func` argument of the `EnsembleProblem` function is used to define the function that generates the problem for each iteration. The `iter` variable contains the product of the `ωq_list` and `ωd_list` lists, which are used to sweep over the parameters. The `sol = solve(ens_prob, Tsit5(), EnsembleDistributed(), trajectories=length(iter))` command is used to solve the problem with the distributed ensemble method. The output of the script will be printed in the `output.out` file, which contains an output similar to the previous example.
+Notice that we are using the [`mesolve_map`](@ref) function, which internally uses the `SciMLBase.EnsembleProblem` function to parallelize the computation of [`mesolveProblem`](@ref). The result is an `Array` of [`TimeEvolutionSol`](@ref), where each element corresponds to a specific combination of initial states and parameters. One can access the solution for a specific combination of initial states and parameters using indexing. For example,
 
+```julia
+sol[i,j,k,l].expect
+```
 
+this will give the expectation values for the case of `ψ_list[i]`, `ωq_list[j]`, `ωd_list[k]`, and `F_list[l]`.