Improve code quality tests

src/QuantumToolbox.jl

@@ -60,7 +60,7 @@ import DiffEqNoiseProcess: RealWienerProcess!
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
 import Distributed: RemoteChannel
-import FFTW: fft, ifft, fftfreq, fftshift
+import FFTW: fft, ifft, fftfreq, fftshift, unsafe_destroy_plan
 import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
 import Pkg

src/correlations.jl

@@ -57,12 +57,12 @@ function correlation_3op_2t(
     kwargs2 = merge((saveat = collect(tlist),), (; kwargs...))
     ρt_list = mesolve(L, ψ0, tlist; kwargs2...).states
 
-    corr = map((t, ρt) -> mesolve(L, C * ρt * A, τlist .+ t, e_ops = [B]; kwargs...).expect[1, :], tlist, ρt_list)
+    corr = Matrix{eltype(ψ0)}(undef, length(tlist), length(τlist))
+    foreach(enumerate(tlist)) do (j, t)
+        return corr[j, :] .= mesolve(L, C * ρt_list[j] * A, τlist .+ t, e_ops = [B]; kwargs...).expect[1, :]
+    end
 
-    # make the output correlation Matrix align with QuTiP
-    # 1st dimension corresponds to tlist
-    # 2nd dimension corresponds to τlist
-    return reduce(vcat, transpose.(corr))
+    return corr
 end
 
 @doc raw"""

src/spectrum.jl

@@ -127,47 +127,47 @@ function _spectrum(
     _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_CType(L), D)) # same as vec(system_identity_matrix)
     _tr_A = transpose(_tr) * spre(A).data
 
-    cache = nothing
     I_cache = I(D^2)
+
+    ω = popfirst!(ωList)
+    cache = init(LinearProblem(L.data - 1im * ω * I_cache, b), solver.alg, kwargs...)
+    sol = solve!(cache)
+    spec[1] = -2 * real(dot(_tr_A, sol.u))
+
     for (idx, ω) in enumerate(ωList)
-        if idx == 1
-            cache = init(LinearProblem(L.data - 1im * ω * I_cache, b), solver.alg, kwargs...)
-            sol = solve!(cache)
-        else
-            cache.A = L.data - 1im * ω * I_cache
-            sol = solve!(cache)
-        end
+        cache.A = L.data - 1im * ω * I_cache
+        sol = solve!(cache)
 
         # trace over the Hilbert space of system (expectation value)
-        spec[idx] = -2 * real(dot(_tr_A, sol.u))
+        spec[idx+1] = -2 * real(dot(_tr_A, sol.u))
     end
 
     return spec
 end
 
 @doc raw"""
-    spectrum_correlation_fft(tlist, corr; inverse=false)
+    spectrum_correlation_fft(tlist, corr; inverse=Val(false))
 
 Calculate the power spectrum corresponding to a two-time correlation function using fast Fourier transform (FFT).
 
 # Parameters
 - `tlist::AbstractVector`: List of times at which the two-time correlation function is given.
 - `corr::AbstractVector`: List of two-time correlations corresponding to the given time point in `tlist`.
-- `inverse::Bool`: Whether to use the inverse Fourier transform or not. Default to `false`.
+- `inverse::Union{Val,Bool}`: Whether to use the inverse Fourier transform or not. Default to `Val(false)`.
 
 # Returns
 - `ωlist`: the list of angular frequencies ``\omega``.
 - `Slist`: the list of the power spectrum corresponding to the angular frequencies in `ωlist`.
 """
-function spectrum_correlation_fft(tlist::AbstractVector, corr::AbstractVector; inverse::Bool = false)
+function spectrum_correlation_fft(tlist::AbstractVector, corr::AbstractVector; inverse::Union{Val,Bool} = Val(false))
     N = length(tlist)
     dt_list = diff(tlist)
     dt = dt_list[1]
 
     all(≈(dt), dt_list) || throw(ArgumentError("tlist must be equally spaced for FFT."))
 
     # power spectrum list
-    F = inverse ? N * ifft(corr) : fft(corr)
+    F = getVal(inverse) ? N * ifft(corr) : fft(corr)
     Slist = 2 * dt * real(fftshift(F))
 
     # angular frequency list

src/steadystate.jl

@@ -167,7 +167,8 @@ function _steadystate(
     Tn = sparse(rows, cols, vals, N^2, N^2)
     L_tmp = L_tmp + Tn
 
-    ρss_vec = L_tmp \ v0 # This is still not supported on GPU, yet
+    fac = lu(L_tmp)
+    ρss_vec = fac \ v0 # This is still not supported on GPU, yet
     ρss = reshape(ρss_vec, N, N)
     ρss = (ρss + ρss') / 2 # Hermitianize
     return QuantumObject(ρss, Operator, L.dims)

src/time_evolution/mesolve.jl

@@ -76,6 +76,7 @@ function mesolveProblem(
     check_dims(L_evo, ψ0)
 
     T = Base.promote_eltype(L_evo, ψ0)
+
     ρ0 = sparse_to_dense(_CType(T), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
     L = L_evo.data
 

test/code-quality/code_quality.jl

@@ -0,0 +1,31 @@
+@testset "Code quality" verbose = true begin
+    @testset "Aqua.jl" begin
+        Aqua.test_all(QuantumToolbox; ambiguities = false, unbound_args = false)
+    end
+
+    @testset "JET.jl" begin
+        # JET.test_package(QuantumToolbox; target_defined_modules = true, ignore_missing_comparison = true)
+
+        include("workload.jl")
+
+        # Here we define some functins to exclude from the JET test
+
+        # Related to the `check_error` function inside the `integrator` interface
+        sci_ml_integrator_functions = (
+            Base.modulesof!,
+            Base.show,
+            Base.show_at_namedtuple,
+            Base.show_typealias,
+            Base._show_type,
+            Base.isvisible,
+            Base.eltype,
+        )
+
+        # Related to FFTW.jl
+        fftw_functions = (QuantumToolbox.unsafe_destroy_plan,)
+
+        function_filter(@nospecialize f) = f ∉ (sci_ml_integrator_functions..., fftw_functions...)
+
+        @test_opt function_filter = function_filter run_all_functions()
+    end
+end

test/code-quality/workload.jl

@@ -0,0 +1,160 @@
+function run_all_functions()
+    # block diagonal form
+    run_block_diagonal_form()
+
+    # correlations and spectrum
+    run_correlations_and_spectrum()
+
+    # dynamical fock dimension
+    # run_dynamical_fock_dimension() # TODO: fix type instabilities here
+
+    # dynamical shifted fock
+    # run_dynamical_shifted_fock() # TODO: fix type instabilities here
+
+    return nothing
+end
+
+#=
+    Block Diagonal Form
+=#
+function run_block_diagonal_form()
+    H, c_ops, a = driven_dissipative_kerr()
+    L = liouvillian(H, c_ops)
+
+    bdf = block_diagonal_form(L)
+
+    return nothing
+end
+
+#=
+    Correlations and Spectrum
+=#
+function run_correlations_and_spectrum()
+    N = 10
+    H, c_ops, a = driven_dissipative_harmonic_oscillator(nth = 0.01, N = N)
+
+    t_l = range(0, 333 * π, length = 1000)
+    corr1 = correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))
+    corr2 = correlation_3op_1t(H, nothing, t_l, c_ops, qeye(a.dims[1]), a', a; progress_bar = Val(false))
+    ω_l1, spec1 = spectrum_correlation_fft(t_l, corr1)
+
+    ω_l2 = range(0, 3, length = 1000)
+    spec2 = spectrum(H, ω_l2, c_ops, a', a)
+    # spec3 = spectrum(H, ω_l2, c_ops, a', a; solver = PseudoInverse())
+
+    return nothing
+end
+
+#=
+    Dynamical Fock Dimension
+=#
+function run_dynamical_fock_dimension()
+    t_l = range(0, 15, length = 100)
+
+    # Single mode
+    F, Δ, κ = 5, 0.25, 1
+
+    maxdims = [150]
+    ψ0 = fock(3, 0)
+    dfd_params = (Δ = Δ, F = F, κ = κ)
+    sol = dfd_mesolve(H_dfd1, ψ0, t_l, c_ops_dfd1, maxdims, dfd_params, e_ops = e_ops_dfd1, progress_bar = Val(false))
+
+    # Two modes
+    F, Δ, κ, J = 1.5, 0.25, 1, 0.05
+
+    maxdims = [50, 50]
+    ψ0 = kron(fock(3, 0), fock(20, 15))
+    dfd_params = (Δ = Δ, F = F, κ = κ, J = J)
+    sol = dfd_mesolve(H_dfd2, ψ0, t_l, c_ops_dfd2, maxdims, dfd_params, e_ops = e_ops_dfd2, progress_bar = Val(false))
+
+    return nothing
+end
+
+#=
+    Dynamical Shifted Fock
+=#
+function run_dynamical_shifted_fock()
+    # Single mode
+    F = 3
+    Δ = 0.25
+    κ = 1
+    U = 0.01
+    α0 = 1.5
+
+    tlist = range(0, 25, 300)
+
+    N = 5
+    a = destroy(N)
+
+    op_list = [a]
+    ψ0 = fock(N, 0)
+    α0_l = [α0]
+    dsf_params = (Δ = Δ, F = F, κ = κ, U = U)
+
+    sol_dsf_me = dsf_mesolve(
+        H_dsf,
+        ψ0,
+        tlist,
+        c_ops_dsf,
+        op_list,
+        α0_l,
+        dsf_params,
+        e_ops = e_ops_dsf,
+        progress_bar = Val(false),
+    )
+    sol_dsf_mc = dsf_mcsolve(
+        H_dsf,
+        ψ0,
+        tlist,
+        c_ops_dsf,
+        op_list,
+        α0_l,
+        dsf_params,
+        e_ops = e_ops_dsf,
+        progress_bar = Val(false),
+        ntraj = 500,
+    )
+
+    # Two modes
+    F = 2
+    Δ = 0.25
+    κ = 1
+    U = 0.01
+    J = 0.5
+    tlist = range(0, 15, 300)
+
+    N = 5
+    a1 = kron(destroy(N), qeye(N))
+    a2 = kron(qeye(N), destroy(N))
+
+    op_list = [a1, a2]
+    ψ0 = kron(fock(N, 0), fock(N, 0))
+    α0_l = [α0, α0]
+    dsf_params = (Δ = Δ, F = F, κ = κ, U = U, J = J)
+
+    sol_dsf_me = dsf_mesolve(
+        H_dsf2,
+        ψ0,
+        tlist,
+        c_ops_dsf2,
+        op_list,
+        α0_l,
+        dsf_params,
+        e_ops = e_ops_dsf2,
+        progress_bar = Val(false),
+    )
+    sol_dsf_mc = dsf_mcsolve(
+        H_dsf2,
+        ψ0,
+        tlist,
+        c_ops_dsf2,
+        op_list,
+        α0_l,
+        dsf_params,
+        e_ops = e_ops_dsf2,
+        progress_bar = Val(false),
+        ntraj = 500,
+    )
+
+    return nothing
+end

test/core-test/block_diagonal_form.jl

@@ -1,19 +1,5 @@
 @testset "Block Diagonal Form" begin
-    # Block Diagonal Form
-    N = 20
-    Δ = 0
-    G = 5
-    tg = 0
-    θ = atan(tg)
-    U = sin(θ)
-    κ2 = cos(θ)
-    κϕ = 1e-3
-    nth = 0.0
-
-    a = destroy(N)
-    ad = create(N)
-    H = -Δ * ad * a + G / 2 * (ad^2 + a^2) + U / 2 * (ad^2 * a^2)
-    c_ops = [√(κ2) * a^2, √(κϕ) * ad * a]
+    H, c_ops, a = driven_dissipative_kerr()
     L = liouvillian(H, c_ops)
 
     bdf = block_diagonal_form(L)

test/core-test/correlations_and_spectrum.jl

@@ -1,9 +1,7 @@
 @testset "Correlations and Spectrum" begin
     N = 10
+    H, c_ops, a = driven_dissipative_harmonic_oscillator(nth = 0.01, N = N)
     Id = qeye(N)
-    a = destroy(N)
-    H = a' * a
-    c_ops = [sqrt(0.1 * (0.01 + 1)) * a, sqrt(0.1 * (0.01)) * a']
 
     t_l = range(0, 333 * π, length = 1000)
     corr1 = correlation_2op_1t(H, nothing, t_l, c_ops, a', a; progress_bar = Val(false))

test/core-test/dynamical-shifted-fock.jl

@@ -18,22 +18,7 @@
 
     N = 5
     a = destroy(N)
-    function H_dsf(op_list, p)
-        Δ = p.Δ
-        F = p.F
-        U = p.U
-        a = op_list[1]
-        return Δ * a' * a + F * (a + a') + U * a'^2 * a^2
-    end
-    function c_ops_dsf(op_list, p)
-        κ = p.κ
-        a = op_list[1]
-        return [√κ * a]
-    end
-    function e_ops_dsf(op_list, p)
-        a = op_list[1]
-        return [a' * a, a]
-    end
+
     op_list = [a]
     ψ0 = fock(N, 0)
     α0_l = [α0]
@@ -92,28 +77,7 @@
     N = 5
     a1 = kron(destroy(N), qeye(N))
     a2 = kron(qeye(N), destroy(N))
-    function H_dsf2(op_list, p)
-        Δ = p.Δ
-        F = p.F
-        U = p.U
-        J = p.J
-        a1, a2 = op_list
-        return Δ * a1' * a1 +
-               Δ * a2' * a2 +
-               U * a1'^2 * a1^2 +
-               U * a2'^2 * a2^2 +
-               F * (a1 + a1') +
-               J * (a1' * a2 + a1 * a2')
-    end
-    function c_ops_dsf2(op_list, p)
-        κ = p.κ
-        a1, a2 = op_list
-        return [√κ * a1, √κ * a2]
-    end
-    function e_ops_dsf2(op_list, p)
-        a1, a2 = op_list
-        return [a1' * a1, a2' * a2]
-    end
+
     op_list = [a1, a2]
     ψ0 = kron(fock(N, 0), fock(N, 0))
     α0_l = [α0, α0]