From 76af7e7ab2d7e5736e2f674ebd7f61157f5293ee Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sat, 5 Apr 2025 17:07:20 +0200 Subject: [PATCH 1/2] Format the code according to the new version of JuliaFormatter.jl --- src/entropy.jl | 2 +- src/qobj/eigsolve.jl | 19 +++++----- src/qobj/operators.jl | 12 +++---- src/qobj/states.jl | 2 +- src/spin_lattice.jl | 6 ++-- src/steadystate.jl | 8 ++--- src/time_evolution/lr_mesolve.jl | 36 +++++++------------ src/time_evolution/time_evolution.jl | 4 +-- .../time_evolution_dynamical.jl | 10 +++--- src/wigner.jl | 6 ++-- test/core-test/generalized_master_equation.jl | 4 +-- test/core-test/low_rank_dynamics.jl | 6 ++-- test/core-test/steady_state.jl | 4 +-- 13 files changed, 53 insertions(+), 66 deletions(-) diff --git a/src/entropy.jl b/src/entropy.jl index 43560fd6d..d32bc7c9a 100644 --- a/src/entropy.jl +++ b/src/entropy.jl @@ -55,7 +55,7 @@ function entropy_vn( length(indexes) == 0 && return zero(real(T)) nzvals = vals[indexes] logvals = base != 0 ? log.(base, Complex.(nzvals)) : log.(Complex.(nzvals)) - return -real(mapreduce(*, +, nzvals, logvals)) + return -real(mapreduce(*,+,nzvals,logvals)) end @doc raw""" diff --git a/src/qobj/eigsolve.jl b/src/qobj/eigsolve.jl index f180c0272..d68a3a163 100644 --- a/src/qobj/eigsolve.jl +++ b/src/qobj/eigsolve.jl @@ -225,7 +225,7 @@ function _eigsolve( # println( A * view(V, :, 1:k) ≈ view(V, :, 1:k) * M(view(H, 1:k, 1:k)) + qₘ * M(transpose(view(transpose(βeₘ) * Uₘ, 1:k))) ) # SHOULD BE TRUE - for j in k+1:m + for j in (k+1):m β = arnoldi_step!(A, V, H, j) if β < tol numops += j - k - 1 @@ -406,15 +406,14 @@ function eigsolve_al( kwargs..., ) where {HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} L_evo = _mesolve_make_L_QobjEvo(H, c_ops) - prob = - mesolveProblem( - L_evo, - QuantumObject(ρ0, type = Operator, dims = H.dimensions), - [zero(T), T]; - params = params, - progress_bar = Val(false), - kwargs..., - ).prob + prob = mesolveProblem( + L_evo, + QuantumObject(ρ0, type = Operator, dims = H.dimensions), + [zero(T), T]; + params = params, + progress_bar = Val(false), + kwargs..., + ).prob integrator = init(prob, alg) Lmap = ArnoldiLindbladIntegratorMap(eltype(H), size(L_evo), integrator) diff --git a/src/qobj/operators.jl b/src/qobj/operators.jl index 5445818ef..90d6d1739 100644 --- a/src/qobj/operators.jl +++ b/src/qobj/operators.jl @@ -100,7 +100,7 @@ julia> fock(20, 3)' * a * fock(20, 4) 2.0 + 0.0im ``` """ -destroy(N::Int) = QuantumObject(spdiagm(1 => Array{ComplexF64}(sqrt.(1:N-1))), Operator, N) +destroy(N::Int) = QuantumObject(spdiagm(1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator, N) @doc raw""" create(N::Int) @@ -126,7 +126,7 @@ julia> fock(20, 4)' * a_d * fock(20, 3) 2.0 + 0.0im ``` """ -create(N::Int) = QuantumObject(spdiagm(-1 => Array{ComplexF64}(sqrt.(1:N-1))), Operator, N) +create(N::Int) = QuantumObject(spdiagm(-1 => Array{ComplexF64}(sqrt.(1:(N-1)))), Operator, N) @doc raw""" displace(N::Int, α::Number) @@ -167,7 +167,7 @@ Bosonic number operator with Hilbert space cutoff `N`. This operator is defined as ``\hat{N}=\hat{a}^\dagger \hat{a}``, where ``\hat{a}`` is the bosonic annihilation operator. """ -num(N::Int) = QuantumObject(spdiagm(0 => Array{ComplexF64}(0:N-1)), Operator, N) +num(N::Int) = QuantumObject(spdiagm(0 => Array{ComplexF64}(0:(N-1))), Operator, N) @doc raw""" position(N::Int) @@ -311,11 +311,11 @@ end jmat(j::Real, ::Val{T}) where {T} = throw(ArgumentError("Invalid spin operator: $(T)")) function _jm(j::Real) - m = j:(-1):-j - data = sqrt.(j * (j + 1) .- m .* (m .- 1))[1:end-1] + m = j:(-1):(-j) + data = sqrt.(j * (j + 1) .- m .* (m .- 1))[1:(end-1)] return spdiagm(-1 => Array{ComplexF64}(data)) end -_jz(j::Real) = spdiagm(0 => Array{ComplexF64}(j .- (0:Int(2 * j)))) +_jz(j::Real) = spdiagm(0 => Array{ComplexF64}(j .- (0:Int(2*j)))) @doc raw""" spin_Jx(j::Real) diff --git a/src/qobj/states.jl b/src/qobj/states.jl index 784c172aa..c2759c9a4 100644 --- a/src/qobj/states.jl +++ b/src/qobj/states.jl @@ -127,7 +127,7 @@ Density matrix for a thermal state (generating thermal state probabilities) with """ function thermal_dm(N::Int, n::Real; sparse::Union{Bool,Val} = Val(false)) β = log(1.0 / n + 1.0) - N_list = Array{Float64}(0:N-1) + N_list = Array{Float64}(0:(N-1)) data = exp.(-β .* N_list) if getVal(sparse) return QuantumObject(spdiagm(0 => data ./ sum(data)), Operator, N) diff --git a/src/spin_lattice.jl b/src/spin_lattice.jl index 17275aba1..72ee4f38f 100644 --- a/src/spin_lattice.jl +++ b/src/spin_lattice.jl @@ -53,11 +53,11 @@ function multisite_operator(dims::Union{AbstractVector,Tuple}, pairs::Pair{<:Int _dims[sites] == [get_dimensions_to(op)[1].size for op in ops] || throw(ArgumentError("The dimensions of the operators do not match the dimensions of the lattice.")) - data = kron(I(prod(_dims[1:sites[1]-1])), ops[1].data) + data = kron(I(prod(_dims[1:(sites[1]-1)])), ops[1].data) for i in 2:length(sites) - data = kron(data, I(prod(_dims[sites[i-1]+1:sites[i]-1])), ops[i].data) + data = kron(data, I(prod(_dims[(sites[i-1]+1):(sites[i]-1)])), ops[i].data) end - data = kron(data, I(prod(_dims[sites[end]+1:end]))) + data = kron(data, I(prod(_dims[(sites[end]+1):end]))) return QuantumObject(data; type = Operator, dims = dims) end diff --git a/src/steadystate.jl b/src/steadystate.jl index 921c7b398..1d25eda6b 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -367,7 +367,7 @@ function _steadystate_fourier( N = size(L_0_mat, 1) Ns = isqrt(N) n_fourier = 2 * n_max + 1 - n_list = -n_max:n_max + n_list = (-n_max):n_max weight = 1 Mn = sparse(ones(Ns), [Ns * (j - 1) + j for j in 1:Ns], fill(weight, Ns), N, N) @@ -399,15 +399,15 @@ function _steadystate_fourier( offset1 = n_max * N offset2 = (n_max + 1) * N - ρ0 = reshape(ρtot[offset1+1:offset2], Ns, Ns) + ρ0 = reshape(ρtot[(offset1+1):offset2], Ns, Ns) ρ0_tr = tr(ρ0) ρ0 = ρ0 / ρ0_tr ρ0 = QuantumObject((ρ0 + ρ0') / 2, type = Operator, dims = L_0.dimensions) ρtot = ρtot / ρ0_tr ρ_list = [ρ0] - for i in 0:n_max-1 - ρi_m = reshape(ρtot[offset1-(i+1)*N+1:offset1-i*N], Ns, Ns) + for i in 0:(n_max-1) + ρi_m = reshape(ρtot[(offset1-(i+1)*N+1):(offset1-i*N)], Ns, Ns) ρi_m = QuantumObject(ρi_m, type = Operator, dims = L_0.dimensions) push!(ρ_list, ρi_m) end diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl index 873e7e9ed..1a06445f6 100644 --- a/src/time_evolution/lr_mesolve.jl +++ b/src/time_evolution/lr_mesolve.jl @@ -55,9 +55,7 @@ lr_mesolve_options_default = ( Δt = 0.0, ) -#=======================================================# # ADDITIONAL FUNCTIONS -#=======================================================# select(x::Real, xarr::AbstractArray, retval = false) = retval ? xarr[argmin(abs.(x .- xarr))] : argmin(abs.(x .- xarr)) @@ -133,9 +131,7 @@ function _calculate_expectation!(p, z, B, idx) end end -#=======================================================# # SAVING FUNCTIONS -#=======================================================# function _periodicsave_func(integrator) ip = integrator.p @@ -151,8 +147,8 @@ function _save_affect_lr_mesolve!(integrator) N, M = ip.N, ip.M idx = select(integrator.t, ip.times) - @views z = reshape(integrator.u[1:N*M], N, M) - @views B = reshape(integrator.u[N*M+1:end], M, M) + @views z = reshape(integrator.u[1:(N*M)], N, M) + @views B = reshape(integrator.u[(N*M+1):end], M, M) _calculate_expectation!(ip, z, B, idx) if integrator.p.opt.progress @@ -162,9 +158,7 @@ function _save_affect_lr_mesolve!(integrator) return u_modified!(integrator, false) end -#=======================================================# # CALLBACK FUNCTIONS -#=======================================================# #= _adjM_condition_ratio(u, t, integrator) @@ -185,8 +179,8 @@ function _adjM_condition_ratio(u, t, integrator) C = ip.A0 σ = ip.temp_MM - @views z = reshape(u[1:N*M], N, M) - @views B = reshape(u[N*M+1:end], M, M) + @views z = reshape(u[1:(N*M)], N, M) + @views B = reshape(u[(N*M+1):end], M, M) mul!(C, z, sqrt(B)) mul!(σ, C', C) p = abs.(eigvals(σ)) @@ -232,8 +226,8 @@ function _adjM_affect!(integrator) N, M = ip.N, ip.M @views begin - z = Δt > 0 ? reshape(ip.u_save[1:N*M], N, M) : reshape(integrator.u[1:N*M], N, M) - B = Δt > 0 ? reshape(ip.u_save[N*M+1:end], M, M) : reshape(integrator.u[N*M+1:end], M, M) + z = Δt > 0 ? reshape(ip.u_save[1:(N*M)], N, M) : reshape(integrator.u[1:(N*M)], N, M) + B = Δt > 0 ? reshape(ip.u_save[(N*M+1):end], M, M) : reshape(integrator.u[(N*M+1):end], M, M) ψ = ip.L_tilde[:, 1] normalize!(ψ) @@ -241,7 +235,7 @@ function _adjM_affect!(integrator) B = cat(B, opt.p0, dims = (1, 2)) resize!(integrator, length(z) + length(B)) integrator.u[1:length(z)] .= z[:] - integrator.u[length(z)+1:end] .= B[:] + integrator.u[(length(z)+1):end] .= B[:] end integrator.p = merge( @@ -277,9 +271,7 @@ function _adjM_affect!(integrator) end end -#=======================================================# # DYNAMICAL EVOLUTION EQUATIONS -#=======================================================# #= dBdz!(du, u, p, t) @@ -305,10 +297,10 @@ function dBdz!(du, u, p, t) N, M = p.N, p.M opt = p.opt - @views z = reshape(u[1:N*M], N, M) - @views dz = reshape(du[1:N*M], N, M) - @views B = reshape(u[N*M+1:end], M, M) - @views dB = reshape(du[N*M+1:end], M, M) + @views z = reshape(u[1:(N*M)], N, M) + @views dz = reshape(du[1:(N*M)], N, M) + @views B = reshape(u[(N*M+1):end], M, M) + @views dB = reshape(du[(N*M+1):end], M, M) #Assign A0 and S mul!(S, z', z) @@ -352,17 +344,13 @@ function dBdz!(du, u, p, t) return dB .-= temp_MM end -#=======================================================# # OUTPUT GENNERATION -#=======================================================# -get_z(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, 1:M*N), N, M) +get_z(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, 1:(M*N)), N, M) get_B(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, (M*N+1):length(u)), M, M) -#=======================================================# # PROBLEM FORMULATION -#=======================================================# @doc raw""" lr_mesolveProblem( diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl index b487d6511..13188cc25 100644 --- a/src/time_evolution/time_evolution.jl +++ b/src/time_evolution/time_evolution.jl @@ -397,7 +397,7 @@ end @generated function update_coefficients!(L::DiffusionOperator, u, p, t) ops_types = L.parameters[2].parameters N = length(ops_types) - quote + return quote Base.@nexprs $N i -> begin update_coefficients!(L.ops[i], u, p, t) end @@ -538,7 +538,7 @@ function _liouvillian_floquet( S = -(L_0 - 1im * n_max * ω * I) \ L_p_dense T = -(L_0 + 1im * n_max * ω * I) \ L_m_dense - for n_i in n_max-1:-1:1 + for n_i in (n_max-1):-1:1 S = -(L_0 - 1im * n_i * ω * I + L_m * S) \ L_p_dense T = -(L_0 + 1im * n_i * ω * I + L_p * T) \ L_m_dense end diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl index c320946e8..e46633239 100644 --- a/src/time_evolution/time_evolution_dynamical.jl +++ b/src/time_evolution/time_evolution_dynamical.jl @@ -39,15 +39,15 @@ function _increase_dims( if n_d == 1 ρmat = similar(QO, dims_new[1], dims_new[1]) - fill!(selectdim(ρmat, 1, dims[1]+1:dims_new[1]), 0) - fill!(selectdim(ρmat, 2, dims[1]+1:dims_new[1]), 0) + fill!(selectdim(ρmat, 1, (dims[1]+1):dims_new[1]), 0) + fill!(selectdim(ρmat, 2, (dims[1]+1):dims_new[1]), 0) copyto!(view(ρmat, 1:dims[1], 1:dims[1]), QO) else ρmat2 = similar(QO, reverse(vcat(dims_new, dims_new))...) ρmat = reshape(QO, reverse(vcat(dims, dims))...) for i in eachindex(sel) - fill!(selectdim(ρmat2, n_d - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0) - fill!(selectdim(ρmat2, 2 * n_d - sel[i] + 1, dims[sel[i]]+1:dims_new[sel[i]]), 0) + fill!(selectdim(ρmat2, n_d - sel[i] + 1, (dims[sel[i]]+1):dims_new[sel[i]]), 0) + fill!(selectdim(ρmat2, 2 * n_d - sel[i] + 1, (dims[sel[i]]+1):dims_new[sel[i]]), 0) end copyto!(view(ρmat2, reverse!(repeat([1:n for n in dims], 2))...), ρmat) ρmat = reshape(ρmat2, prod(dims_new), prod(dims_new)) @@ -77,7 +77,7 @@ function _DFDIncreaseReduceCondition(u, t, integrator) pillow_i = pillow_list[i] if dim_i < maxdim_i && dim_i > 2 && maxdim_i != 0 ρi = _ptrace_oper(vec2mat(dfd_ρt_cache), dim_list, SVector(i))[1] - @views res = norm(ρi[diagind(ρi)[end-pillow_i:end]], 1) * sqrt(dim_i) / pillow_i + @views res = norm(ρi[diagind(ρi)[(end-pillow_i):end]], 1) * sqrt(dim_i) / pillow_i if res > tol_list[i] increase_list[i] = true elseif res < tol_list[i] * 1e-2 && dim_i > 3 diff --git a/src/wigner.jl b/src/wigner.jl index dec636642..813457b6e 100644 --- a/src/wigner.jl +++ b/src/wigner.jl @@ -166,10 +166,10 @@ function _wigner_laguerre(ρ::AbstractArray, A::AbstractArray, W::AbstractArray, if method.parallel throw(ArgumentError("Parallel version is not implemented for dense matrices")) else - for m in 0:M-1 + for m in 0:(M-1) ρmn = ρ[m+1, m+1] abs(ρmn) > tol && (@. W += real(ρmn * (-1)^m * _genlaguerre(m, 0, B))) - for n in m+1:M-1 + for n in (m+1):(M-1) ρmn = ρ[m+1, n+1] # Γ_mn = sqrt(gamma(m+1) / gamma(n+1)) Γ_mn = sqrt(exp(loggamma(m + 1) - loggamma(n + 1))) # Is this a good trick? @@ -197,7 +197,7 @@ function _genlaguerre(n::Int, α::Number, x::T) where {T<:BlasFloat} α = convert(T, α) p0, p1 = one(T), -x + (α + 1) n == 0 && return p0 - for k in 1:n-1 + for k in 1:(n-1) p1, p0 = ((2k + α + 1) / (k + 1) - x / (k + 1)) * p1 - (k + α) / (k + 1) * p0, p1 end return p1 diff --git a/test/core-test/generalized_master_equation.jl b/test/core-test/generalized_master_equation.jl index 54a9ac8e9..03d2989e0 100644 --- a/test/core-test/generalized_master_equation.jl +++ b/test/core-test/generalized_master_equation.jl @@ -15,7 +15,7 @@ Tlist = [0, 0.0] E, U, L1 = liouvillian_generalized(H, fields, Tlist, N_trunc = N_trunc, tol = tol) - Ω = to_sparse((E'.-E)[1:N_trunc, 1:N_trunc], tol) + Ω = to_sparse((E' .- E)[1:N_trunc, 1:N_trunc], tol) H_d = Qobj(to_sparse((U'*H*U)[1:N_trunc, 1:N_trunc], tol)) Xp = Qobj(Ω .* to_sparse(triu((U'*(a+a')*U).data[1:N_trunc, 1:N_trunc], 1), tol)) @@ -33,7 +33,7 @@ Tlist = [0.2, 0.0] E, U, L1 = liouvillian_generalized(H, fields, Tlist, N_trunc = N_trunc, tol = tol) - Ω = to_sparse((E'.-E)[1:N_trunc, 1:N_trunc], tol) + Ω = to_sparse((E' .- E)[1:N_trunc, 1:N_trunc], tol) H_d = Qobj(to_sparse((U'*H*U)[1:N_trunc, 1:N_trunc], tol)) Xp = Qobj(Ω .* to_sparse(triu((U'*(a+a')*U).data[1:N_trunc, 1:N_trunc], 1), tol)) diff --git a/test/core-test/low_rank_dynamics.jl b/test/core-test/low_rank_dynamics.jl index 7b718c734..27066cdb1 100644 --- a/test/core-test/low_rank_dynamics.jl +++ b/test/core-test/low_rank_dynamics.jl @@ -17,13 +17,13 @@ i += 1 i <= M && (ϕ[i] = multisite_operator(latt, j => sigmap()) * ϕ[1]) end - for k in 1:N_modes-1 - for l in k+1:N_modes + for k in 1:(N_modes-1) + for l in (k+1):N_modes i += 1 i <= M && (ϕ[i] = multisite_operator(latt, k => sigmap(), l => sigmap()) * ϕ[1]) end end - for i in i+1:M + for i in (i+1):M ϕ[i] = QuantumObject(rand(ComplexF64, size(ϕ[1])[1]), dims = ϕ[1].dims) normalize!(ϕ[i]) end diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl index 795ef7d5a..9981d5abb 100644 --- a/test/core-test/steady_state.jl +++ b/test/core-test/steady_state.jl @@ -68,8 +68,8 @@ ρ_ss2 = steadystate_fourier(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian()) - @test abs(sum(sol_me.expect[1, end-100:end]) / 101 - expect(e_ops[1], ρ_ss1)) < 1e-3 - @test abs(sum(sol_me.expect[1, end-100:end]) / 101 - expect(e_ops[1], ρ_ss2)) < 1e-3 + @test abs(sum(sol_me.expect[1, (end-100):end]) / 101 - expect(e_ops[1], ρ_ss1)) < 1e-3 + @test abs(sum(sol_me.expect[1, (end-100):end]) / 101 - expect(e_ops[1], ρ_ss2)) < 1e-3 @testset "Type Inference (steadystate_fourier)" begin @inferred steadystate_fourier( From b79f8b6b3e70186168ab1a6777823b2788695c75 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Sun, 6 Apr 2025 10:53:53 +0200 Subject: [PATCH 2/2] Make multiline comments --- src/time_evolution/lr_mesolve.jl | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl index 1a06445f6..4a0f204b3 100644 --- a/src/time_evolution/lr_mesolve.jl +++ b/src/time_evolution/lr_mesolve.jl @@ -55,7 +55,9 @@ lr_mesolve_options_default = ( Δt = 0.0, ) -# ADDITIONAL FUNCTIONS +#= + ADDITIONAL FUNCTIONS +=# select(x::Real, xarr::AbstractArray, retval = false) = retval ? xarr[argmin(abs.(x .- xarr))] : argmin(abs.(x .- xarr)) @@ -131,7 +133,9 @@ function _calculate_expectation!(p, z, B, idx) end end -# SAVING FUNCTIONS +#= + SAVING FUNCTIONS +=# function _periodicsave_func(integrator) ip = integrator.p @@ -158,7 +162,9 @@ function _save_affect_lr_mesolve!(integrator) return u_modified!(integrator, false) end -# CALLBACK FUNCTIONS +#= + CALLBACK FUNCTIONS +=# #= _adjM_condition_ratio(u, t, integrator) @@ -271,7 +277,9 @@ function _adjM_affect!(integrator) end end -# DYNAMICAL EVOLUTION EQUATIONS +#= + DYNAMICAL EVOLUTION EQUATIONS +=# #= dBdz!(du, u, p, t) @@ -344,13 +352,17 @@ function dBdz!(du, u, p, t) return dB .-= temp_MM end -# OUTPUT GENNERATION +#= + OUTPUT GENNERATION +=# get_z(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, 1:(M*N)), N, M) get_B(u::AbstractArray{T}, N::Integer, M::Integer) where {T} = reshape(view(u, (M*N+1):length(u)), M, M) -# PROBLEM FORMULATION +#= + PROBLEM FORMULATION +=# @doc raw""" lr_mesolveProblem(