From 5de07e8bdffdb670d945ec59db716d7ca8462cd6 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sun, 22 Dec 2024 15:50:39 +0100 Subject: [PATCH] Improve code quality tests --- src/QuantumToolbox.jl | 2 +- src/correlations.jl | 10 +- src/spectrum.jl | 26 +-- src/steadystate.jl | 3 +- src/time_evolution/mesolve.jl | 1 + test/code-quality/code_quality.jl | 31 ++++ test/code-quality/workload.jl | 160 ++++++++++++++++++ test/core-test/block_diagonal_form.jl | 16 +- test/core-test/code_quality.jl | 9 - test/core-test/correlations_and_spectrum.jl | 4 +- test/core-test/dynamical-shifted-fock.jl | 40 +---- .../dynamical_fock_dimension_mesolve.jl | 58 +------ test/helpers.jl | 105 ++++++++++++ test/runtests.jl | 5 +- 14 files changed, 330 insertions(+), 140 deletions(-) create mode 100644 test/code-quality/code_quality.jl create mode 100644 test/code-quality/workload.jl delete mode 100644 test/core-test/code_quality.jl create mode 100644 test/helpers.jl diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index 28a50f25..8329f192 100644 --- a/src/QuantumToolbox.jl +++ b/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 diff --git a/src/correlations.jl b/src/correlations.jl index 5b7068d9..800bd707 100644 --- a/src/correlations.jl +++ b/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""" diff --git a/src/spectrum.jl b/src/spectrum.jl index a62a69ca..f486edc4 100644 --- a/src/spectrum.jl +++ b/src/spectrum.jl @@ -127,39 +127,39 @@ 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] @@ -167,7 +167,7 @@ function spectrum_correlation_fft(tlist::AbstractVector, corr::AbstractVector; i 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 diff --git a/src/steadystate.jl b/src/steadystate.jl index 606b0907..4f5b8a6a 100644 --- a/src/steadystate.jl +++ b/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) diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl index 2d9f42ef..99258f35 100644 --- a/src/time_evolution/mesolve.jl +++ b/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 diff --git a/test/code-quality/code_quality.jl b/test/code-quality/code_quality.jl new file mode 100644 index 00000000..1b78c833 --- /dev/null +++ b/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 diff --git a/test/code-quality/workload.jl b/test/code-quality/workload.jl new file mode 100644 index 00000000..2144eb2f --- /dev/null +++ b/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 diff --git a/test/core-test/block_diagonal_form.jl b/test/core-test/block_diagonal_form.jl index cc96a48d..95d85a96 100644 --- a/test/core-test/block_diagonal_form.jl +++ b/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) diff --git a/test/core-test/code_quality.jl b/test/core-test/code_quality.jl deleted file mode 100644 index 5effb97d..00000000 --- a/test/core-test/code_quality.jl +++ /dev/null @@ -1,9 +0,0 @@ -@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) - end -end diff --git a/test/core-test/correlations_and_spectrum.jl b/test/core-test/correlations_and_spectrum.jl index 381cd9e0..d3d2e159 100644 --- a/test/core-test/correlations_and_spectrum.jl +++ b/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)) diff --git a/test/core-test/dynamical-shifted-fock.jl b/test/core-test/dynamical-shifted-fock.jl index bb1bebc1..fa57c875 100644 --- a/test/core-test/dynamical-shifted-fock.jl +++ b/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] diff --git a/test/core-test/dynamical_fock_dimension_mesolve.jl b/test/core-test/dynamical_fock_dimension_mesolve.jl index 9c6c8a02..a3af906b 100644 --- a/test/core-test/dynamical_fock_dimension_mesolve.jl +++ b/test/core-test/dynamical_fock_dimension_mesolve.jl @@ -1,35 +1,19 @@ ### DYNAMICAL FOCK DIMENSION ### @testset "Dynamical Fock Dimension" begin F, Δ, κ = 5, 0.25, 1 + N0 = 140 t_l = range(0, 15, length = 100) - N0 = 140 - a0 = destroy(N0) - H0 = Δ * a0' * a0 + F * (a0 + a0') - c_ops0 = [√κ * a0] + H0, c_ops0, a0 = driven_dissipative_harmonic_oscillator(Δ = Δ, γ = κ, F = F, N = N0) e_ops0 = [a0' * a0] + ψ00 = fock(N0, 0) sol0 = mesolve(H0, ψ00, t_l, c_ops0, e_ops = e_ops0, progress_bar = Val(false)) - function H_dfd0(dims, p) - Δ = p.Δ - F = p.F - a = destroy(dims[1]) - return Δ * a' * a + F * (a + a') - end - function c_ops_dfd0(dims, p) - κ = p.κ - a = destroy(dims[1]) - return [√κ * a] - end - function e_ops_dfd0(dims, p) - a = destroy(dims[1]) - return [a' * a] - end maxdims = [150] ψ0 = fock(3, 0) dfd_params = (Δ = Δ, F = F, κ = κ) - sol = dfd_mesolve(H_dfd0, ψ0, t_l, c_ops_dfd0, maxdims, dfd_params, e_ops = e_ops_dfd0, progress_bar = Val(false)) + sol = dfd_mesolve(H_dfd1, ψ0, t_l, c_ops_dfd1, maxdims, dfd_params, e_ops = e_ops_dfd1, progress_bar = Val(false)) @test sum(abs.((sol.expect[1, :] .- sol0.expect[1, :]) ./ (sol0.expect[1, :] .+ 1e-16))) < 0.01 @@ -42,21 +26,6 @@ ψ00 = fock(N0, 50) sol0 = mesolve(H0, ψ00, t_l, c_ops0, e_ops = e_ops0, progress_bar = Val(false)) - function H_dfd1(dims, p) - Δ = p.Δ - F = p.F - a = destroy(dims[1]) - return Δ * a' * a + F * (a + a') - end - function c_ops_dfd1(dims, p) - κ = p.κ - a = destroy(dims[1]) - return [√κ * a] - end - function e_ops_dfd1(dims, p) - a = destroy(dims[1]) - return [a' * a] - end maxdims = [150] ψ0 = fock(70, 50) dfd_params = (Δ = Δ, F = F, κ = κ) @@ -77,25 +46,6 @@ ψ00 = kron(fock(N0, 0), fock(N1, 15)) sol0 = mesolve(H0, ψ00, t_l, c_ops0, e_ops = e_ops0, progress_bar = Val(false)) - function H_dfd2(dims, p) - Δ = p.Δ - F = p.F - J = p.J - a = kron(destroy(dims[1]), qeye(dims[2])) - b = kron(qeye(dims[1]), destroy(dims[2])) - return Δ * a' * a + F * (a + a') + Δ * b' * b + J * (a' * b + a * b') - end - function c_ops_dfd2(dims, p) - κ = p.κ - a = kron(destroy(dims[1]), qeye(dims[2])) - b = kron(qeye(dims[1]), destroy(dims[2])) - return [√κ * a, √κ * b] - end - function e_ops_dfd2(dims, p) - a = kron(destroy(dims[1]), qeye(dims[2])) - b = kron(qeye(dims[1]), destroy(dims[2])) - return [a' * a, b' * b] - end maxdims = [50, 50] ψ0 = kron(fock(3, 0), fock(20, 15)) dfd_params = (Δ = Δ, F = F, κ = κ, J = J) diff --git a/test/helpers.jl b/test/helpers.jl new file mode 100644 index 00000000..f4b07ef2 --- /dev/null +++ b/test/helpers.jl @@ -0,0 +1,105 @@ +function driven_dissipative_harmonic_oscillator(; Δ = 1, γ = 0.1, nth = 0, F = 0, N = 10) + a = destroy(N) + H = Δ * a' * a + F * (a + a') + c_ops = (sqrt(γ * (nth + 1)) * a, sqrt(γ * nth) * a') + + return H, c_ops, a +end + +function driven_dissipative_kerr() + N = 20 + Δ = 0 + G = 5 + tg = 0 + θ = atan(tg) + U = sin(θ) + κ2 = cos(θ) + κϕ = 1e-3 + + 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) + + return H, c_ops, a +end + +#= + Dynamical Fock Dimension +=# +function H_dfd1(dims, p) + Δ = p.Δ + F = p.F + a = destroy(dims[1]) + return Δ * a' * a + F * (a + a') +end +function c_ops_dfd1(dims, p) + κ = p.κ + a = destroy(dims[1]) + return [√κ * a] +end +function e_ops_dfd1(dims, p) + a = destroy(dims[1]) + return [a' * a] +end +function H_dfd2(dims, p) + Δ = p.Δ + F = p.F + J = p.J + a = kron(destroy(dims[1]), qeye(dims[2])) + b = kron(qeye(dims[1]), destroy(dims[2])) + return Δ * a' * a + F * (a + a') + Δ * b' * b + J * (a' * b + a * b') +end +function c_ops_dfd2(dims, p) + κ = p.κ + a = kron(destroy(dims[1]), qeye(dims[2])) + b = kron(qeye(dims[1]), destroy(dims[2])) + return [√κ * a, √κ * b] +end +function e_ops_dfd2(dims, p) + a = kron(destroy(dims[1]), qeye(dims[2])) + b = kron(qeye(dims[1]), destroy(dims[2])) + return [a' * a, b' * b] +end + +#= + Dynamical Shifted Fock +=# +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 +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 diff --git a/test/runtests.jl b/test/runtests.jl index 9081ad84..56015935 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,11 +26,14 @@ core_tests = [ "wigner.jl", ] +# Include helper functions to avoid repetition +include("helpers.jl") + if (GROUP == "All") || (GROUP == "Code-Quality") using QuantumToolbox using Aqua, JET - include(joinpath(testdir, "core-test", "code_quality.jl")) + include(joinpath(testdir, "code-quality", "code_quality.jl")) end if (GROUP == "All") || (GROUP == "Core")