Other fixes and format

Author: albertomercurio

src/linear_maps.jl

@@ -48,7 +48,7 @@ where `integrator` is the ODE integrator for the time-evolution. In this way, we
 """
 abstract type AbstractLinearMap{T,TS} end
 
-Base.eltype(A::AbstractLinearMap{T}) where T = T
+Base.eltype(A::AbstractLinearMap{T}) where {T} = T
 
 Base.size(A::AbstractLinearMap) = A.size
 Base.size(A::AbstractLinearMap, i::Int) = A.size[i]

src/metrics.jl

@@ -53,12 +53,12 @@ function entropy_vn(
     base::Int = 0,
     tol::Real = 1e-15,
 ) where {T}
-    vals = eigvals(ρ)
-    indexes = abs.(vals) .> tol
-    1 ∉ indexes && return 0
+    vals = eigenenergies(ρ)
+    indexes = findall(x -> abs(x) > tol, vals)
+    length(indexes) == 0 && return zero(real(T))
     nzvals = vals[indexes]
     logvals = base != 0 ? log.(base, Complex.(nzvals)) : log.(Complex.(nzvals))
-    return -real(sum(nzvals .* logvals))
+    return -real(mapreduce(*, +, nzvals, logvals))
 end
 
 """

src/qobj/arithmetic_and_attributes.jl

@@ -240,7 +240,9 @@ julia> tr(a' * a)
 LinearAlgebra.tr(
     A::QuantumObject{<:AbstractArray{T},OpType},
 ) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}} = tr(A.data)
-LinearAlgebra.tr(A::QuantumObject{<:Union{<:Hermitian{TF}, Symmetric{TR}},OpType}) where {TF<:BlasFloat,TR<:Real,OpType<:OperatorQuantumObject} = real(tr(A.data))
+LinearAlgebra.tr(
+    A::QuantumObject{<:Union{<:Hermitian{TF},Symmetric{TR}},OpType},
+) where {TF<:BlasFloat,TR<:Real,OpType<:OperatorQuantumObject} = real(tr(A.data))
 
 @doc raw"""
     svdvals(A::QuantumObject)
@@ -515,7 +517,7 @@ function ptrace(QO::QuantumObject{<:AbstractArray,KetQuantumObject}, sel::Union{
     length(QO.dims) == 1 && return QO
 
     ρtr, dkeep = _ptrace_ket(QO.data, QO.dims, SVector(sel))
-    return QuantumObject(ρtr, type=Operator, dims = dkeep)
+    return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 
 ptrace(QO::QuantumObject{<:AbstractArray,BraQuantumObject}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
@@ -524,7 +526,7 @@ function ptrace(QO::QuantumObject{<:AbstractArray,OperatorQuantumObject}, sel::U
     length(QO.dims) == 1 && return QO
 
     ρtr, dkeep = _ptrace_oper(QO.data, QO.dims, SVector(sel))
-    return QuantumObject(ρtr, type=Operator, dims = dkeep)
+    return QuantumObject(ρtr, type = Operator, dims = dkeep)
 end
 ptrace(QO::QuantumObject, sel::Int) = ptrace(QO, SVector(sel))
 

src/qobj/eigsolve.jl

@@ -389,7 +389,7 @@ function eigsolve_al(
 
     # prog = ProgressUnknown(desc="Applications:", showspeed = true, enabled=progress)
 
-    Lmap = ArnoldiLindbladIntegratorMap(eltype(MT1),size(L),integrator)
+    Lmap = ArnoldiLindbladIntegratorMap(eltype(MT1), size(L), integrator)
 
     res = _eigsolve(Lmap, mat2vec(ρ0), L.type, L.dims, k, krylovdim, maxiter = maxiter, tol = eigstol)
     # finish!(prog)

src/qobj/operators.jl

@@ -31,7 +31,8 @@ The `distribution` specifies which of the method used to obtain the unitary matr
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `rand_unitary(dimensions, Val(distribution))` instead of `rand_unitary(dimensions, distribution)`. Also, put `dimensions` as `Tuple` or `SVector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-rand_unitary(dimensions::Int, distribution::Union{Symbol,Val} = Val(:haar)) = rand_unitary(SVector(dimensions), makeVal(distribution))
+rand_unitary(dimensions::Int, distribution::Union{Symbol,Val} = Val(:haar)) =
+    rand_unitary(SVector(dimensions), makeVal(distribution))
 rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, distribution::Union{Symbol,Val} = Val(:haar)) =
     rand_unitary(dimensions, makeVal(distribution))
 function rand_unitary(dimensions::Union{AbstractVector{Int},Tuple}, ::Val{:haar})
@@ -484,7 +485,7 @@ end
 
 Generates the projection operator ``\hat{O} = \dyad{i}{j}`` with Hilbert space dimension `N`.
 """
-projection(N::Int, i::Int, j::Int) = QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type=Operator)
+projection(N::Int, i::Int, j::Int) = QuantumObject(sparse([i + 1], [j + 1], [1.0 + 0.0im], N, N), type = Operator)
 
 @doc raw"""
     tunneling(N::Int, m::Int=1; sparse::Union{Bool,Val{<:Bool}}=Val(false))

src/qobj/states.jl

@@ -33,12 +33,7 @@ It is also possible to specify the list of dimensions `dims` if different subsys
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `fock(N, j, dims=dims, sparse=Val(sparse))` instead of `fock(N, j, dims=dims, sparse=sparse)`. Consider also to use `dims` as a `Tuple` or `SVector` instead of `Vector`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-function fock(
-    N::Int,
-    j::Int = 0;
-    dims::Union{Int,AbstractVector{Int},Tuple} = N,
-    sparse::Union{Bool,Val} = Val(false),
-)
+function fock(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple} = N, sparse::Union{Bool,Val} = Val(false))
     if getVal(makeVal(sparse))
         array = sparsevec([j + 1], [1.0 + 0im], N)
     else

src/steadystate.jl

@@ -31,10 +31,14 @@ end
 
 struct SteadyStateDirectSolver <: SteadyStateSolver end
 struct SteadyStateEigenSolver <: SteadyStateSolver end
-Base.@kwdef struct SteadyStateLinearSolver{MT<:Union{SciMLLinearSolveAlgorithm,Nothing}} <: SteadyStateSolver
-    alg::MT = nothing
-    Pl::Union{Function,Nothing} = nothing
-    Pr::Union{Function,Nothing} = nothing
+Base.@kwdef struct SteadyStateLinearSolver{
+    MT<:Union{SciMLLinearSolveAlgorithm,Nothing},
+    PlT<:Union{Function,Nothing},
+    PrT<:Union{Function,Nothing},
+} <: SteadyStateSolver
+    alg::MT = KrylovJL_GMRES()
+    Pl::PlT = nothing
+    Pr::PrT = nothing
 end
 
 @doc raw"""
@@ -363,13 +367,13 @@ function _steadystate_floquet(
     ρ0 = reshape(ρtot[offset1+1:offset2], Ns, Ns)
     ρ0_tr = tr(ρ0)
     ρ0 = ρ0 / ρ0_tr
-    ρ0 = QuantumObject((ρ0 + ρ0') / 2, dims = L_0.dims)
+    ρ0 = QuantumObject((ρ0 + ρ0') / 2, type = Operator, dims = L_0.dims)
     ρ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)
-        ρi_m = QuantumObject(ρi_m, dims = L_0.dims)
+        ρi_m = QuantumObject(ρi_m, type = Operator, dims = L_0.dims)
         push!(ρ_list, ρi_m)
     end
 

src/time_evolution/time_evolution.jl

@@ -197,10 +197,10 @@ end
 
 function liouvillian_floquet(
     H::QuantumObject{<:AbstractArray{T1},OpType1},
-    c_ops::AbstractVector,
     Hₚ::QuantumObject{<:AbstractArray{T2},OpType2},
     Hₘ::QuantumObject{<:AbstractArray{T3},OpType3},
-    ω::Real;
+    ω::Real,
+    c_ops::Union{AbstractVector,Nothing} = nothing;
     n_max::Int = 3,
     tol::Real = 1e-15,
 ) where {
@@ -226,30 +226,31 @@ function liouvillian_generalized(
     H::QuantumObject{MT,OperatorQuantumObject},
     fields::Vector,
     T_list::Vector{<:Real};
-    N_trunc::Int = size(H, 1),
+    N_trunc::Union{Int,Nothing} = nothing,
     tol::Real = 1e-12,
     σ_filter::Union{Nothing,Real} = nothing,
 ) where {MT<:AbstractMatrix}
     (length(fields) == length(T_list)) || throw(DimensionMismatch("The number of fields, ωs and Ts must be the same."))
 
-    dims = N_trunc == size(H, 1) ? H.dims : SVector(N_trunc)
+    dims = (N_trunc isa Nothing) ? H.dims : SVector(N_trunc)
+    final_size = prod(dims)
     result = eigen(H)
-    E = real.(result.values[1:N_trunc])
+    E = real.(result.values[1:final_size])
     U = QuantumObject(result.vectors, result.type, result.dims)
 
-    H_d = QuantumObject(Diagonal(complex(E)), dims = dims)
+    H_d = QuantumObject(Diagonal(complex(E)), type = Operator, dims = dims)
 
     Ω = E' .- E
     Ωp = triu(dense_to_sparse(Ω, tol), 1)
 
     # Filter in the Hilbert space
     σ = isnothing(σ_filter) ? 500 * maximum([norm(field) / length(field) for field in fields]) : σ_filter
-    F1 = QuantumObject(gaussian.(Ω, 0, σ), dims = dims)
+    F1 = QuantumObject(gaussian.(Ω, 0, σ), type = Operator, dims = dims)
     F1 = dense_to_sparse(F1, tol)
 
     # Filter in the Liouville space
-    # M1 = ones(N_trunc, N_trunc)
-    M1 = similar(E, N_trunc, N_trunc)
+    # M1 = ones(final_size, final_size)
+    M1 = similar(E, final_size, final_size)
     M1 .= 1
     Ω1 = kron(Ω, M1)
     Ω2 = kron(M1, Ω)
@@ -261,16 +262,16 @@ function liouvillian_generalized(
 
     for i in eachindex(fields)
         # The operator that couples the system to the bath in the eigenbasis
-        X_op = dense_to_sparse((U'*fields[i]*U).data[1:N_trunc, 1:N_trunc], tol)
+        X_op = dense_to_sparse((U'*fields[i]*U).data[1:final_size, 1:final_size], tol)
         if ishermitian(fields[i])
             X_op = (X_op + X_op') / 2 # Make sure it's hermitian
         end
 
         # Ohmic reservoir
         N_th = n_th.(Ωp, T_list[i])
-        Sp₀ = QuantumObject(triu(X_op, 1), dims = dims)
-        Sp₁ = QuantumObject(droptol!((@. Ωp * N_th * Sp₀.data), tol), dims = dims)
-        Sp₂ = QuantumObject(droptol!((@. Ωp * (1 + N_th) * Sp₀.data), tol), dims = dims)
+        Sp₀ = QuantumObject(triu(X_op, 1), type = Operator, dims = dims)
+        Sp₁ = QuantumObject(droptol!((@. Ωp * N_th * Sp₀.data), tol), type = Operator, dims = dims)
+        Sp₂ = QuantumObject(droptol!((@. Ωp * (1 + N_th) * Sp₀.data), tol), type = Operator, dims = dims)
         # S0 = QuantumObject( spdiagm(diag(X_op)), dims=dims )
 
         L +=

test/eigenvalues_and_operators.jl

@@ -93,25 +93,25 @@
     @test isapprox(vec2mat(vecs[1]).data * exp(-1im * angle(vecs[1][1])), vec2mat(state3[1]).data, atol = 1e-5)
 
     @testset "Type Inference (eigen)" begin
-      N = 5
-      a = kron(destroy(N), qeye(N))
-      a_d = a'
-      b = kron(qeye(N), destroy(N))
-      b_d = b'
-  
-      ωc = 1
-      ωb = 1
-      g = 0.01
-      κ = 0.1
-      n_thermal = 0.01
-  
-      H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
-      c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
-      L = liouvillian(H, c_ops)
-
-      @inferred eigenstates(H, sparse = false)
-      @inferred eigenstates(H, sparse = true)
-      @inferred eigenstates(L, sparse = true)
-      @inferred eigsolve_al(L, 1 \ (40 * κ), k = 10)
+        N = 5
+        a = kron(destroy(N), qeye(N))
+        a_d = a'
+        b = kron(qeye(N), destroy(N))
+        b_d = b'
+
+        ωc = 1
+        ωb = 1
+        g = 0.01
+        κ = 0.1
+        n_thermal = 0.01
+
+        H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
+        c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
+        L = liouvillian(H, c_ops)
+
+        @inferred eigenstates(H, sparse = false)
+        @inferred eigenstates(H, sparse = true)
+        @inferred eigenstates(L, sparse = true)
+        @inferred eigsolve_al(L, 1 \ (40 * κ), k = 10)
     end
 end

test/entanglement.jl

@@ -5,4 +5,9 @@
     rho = state * state'
     @test entanglement(state, 1) / log(2) ≈ 1
     @test entanglement(rho, 1) / log(2) ≈ 1
+
+    @testset "Type Stability (entanglement)" begin
+        @inferred entanglement(state, 1)
+        @inferred entanglement(rho, 1)
+    end
 end