Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
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
36 changes: 12 additions & 24 deletions src/time_evolution/lr_mesolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,7 @@
Δt = 0.0,
)

#=======================================================#
# 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 +131,7 @@
end
end

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

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

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

#=
_adjM_condition_ratio(u, t, integrator)
Expand All @@ -185,8 +179,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 183 in src/time_evolution/lr_mesolve.jl

View check run for this annotation

Codecov / codecov/patch

src/time_evolution/lr_mesolve.jl#L182-L183

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

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

#=
dBdz!(du, u, p, t)
Expand All @@ -305,10 +297,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 +344,13 @@
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(
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