Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/entropy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"""
Expand Down
19 changes: 9 additions & 10 deletions src/qobj/eigsolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
12 changes: 6 additions & 6 deletions src/qobj/operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/qobj/states.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions src/spin_lattice.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions src/steadystate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
60 changes: 30 additions & 30 deletions src/time_evolution/lr_mesolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,9 @@
Δt = 0.0,
)

#=======================================================#
# ADDITIONAL FUNCTIONS
#=======================================================#
#=
ADDITIONAL FUNCTIONS
=#

select(x::Real, xarr::AbstractArray, retval = false) = retval ? xarr[argmin(abs.(x .- xarr))] : argmin(abs.(x .- xarr))

Expand Down Expand Up @@ -133,9 +133,9 @@
end
end

#=======================================================#
# SAVING FUNCTIONS
#=======================================================#
#=
SAVING FUNCTIONS
=#

function _periodicsave_func(integrator)
ip = integrator.p
Expand All @@ -151,8 +151,8 @@
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
Expand All @@ -162,9 +162,9 @@
return u_modified!(integrator, false)
end

#=======================================================#
# CALLBACK FUNCTIONS
#=======================================================#
#=
CALLBACK FUNCTIONS
=#

#=
_adjM_condition_ratio(u, t, integrator)
Expand All @@ -185,8 +185,8 @@

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)

Check warning on line 189 in src/time_evolution/lr_mesolve.jl

View check run for this annotation

Codecov / codecov/patch

src/time_evolution/lr_mesolve.jl#L188-L189

Added lines #L188 - L189 were not covered by tests
mul!(C, z, sqrt(B))
mul!(σ, C', C)
p = abs.(eigvals(σ))
Expand Down Expand Up @@ -232,16 +232,16 @@
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!(ψ)

z = hcat(z, ψ)
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(
Expand Down Expand Up @@ -277,9 +277,9 @@
end
end

#=======================================================#
# DYNAMICAL EVOLUTION EQUATIONS
#=======================================================#
#=
DYNAMICAL EVOLUTION EQUATIONS
=#

#=
dBdz!(du, u, p, t)
Expand All @@ -305,10 +305,10 @@
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)
Expand Down Expand Up @@ -352,17 +352,17 @@
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_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(
Expand Down
4 changes: 2 additions & 2 deletions src/time_evolution/time_evolution.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions src/time_evolution/time_evolution_dynamical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/wigner.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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?
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions test/core-test/generalized_master_equation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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))
Expand Down
6 changes: 3 additions & 3 deletions test/core-test/low_rank_dynamics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions test/core-test/steady_state.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
Loading