Skip to content

Commit

Permalink
Added OperatorSum into Dynamical Shifted Fock algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@albertomercurio
albertomercurio committed Feb 20, 2024
1 parent 684ddb7 commit 3c53b1a
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 60 deletions.
39 changes: 4 additions & 35 deletions src/time_evolution/time_evolution.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,41 +53,7 @@ end

ContinuousLindbladJumpCallback(;interp_points::Int=10) = ContinuousLindbladJumpCallback(interp_points)

## Time dependent sum of operators

# mutable struct TimeDependentOperatorSum{MT,CT,CFT}
# operators::MT
# coefficients::CT
# coefficient_functions::CFT
# end

# function TimeDependentOperatorSum(coefficient_functions, operators::Vector; params=nothing, init_time=0.0)
# # promote the type of the coefficients and the operators. Remember that the coefficient_functions si a vector of functions and the operators is a vector of QuantumObjects
# T = promote_type(mapreduce(x->eltype(x.data), promote_type, operators),
# mapreduce(eltype, promote_type, [f(init_time,params) for f in coefficient_functions]))

# coefficients = T[f(init_time, params) for f in coefficient_functions]
# return TimeDependentOperatorSum(operators, coefficients, coefficient_functions)
# end

# function update_coefficients!(A::TimeDependentOperatorSum, t, params)
# @inbounds @simd for i in 1:length(A.coefficient_functions)
# A.coefficients[i] = A.coefficient_functions[i](t, params)
# end
# end

# (A::TimeDependentOperatorSum)(t, params) = (update_coefficients!(A, t, params); A)

# @inline function LinearAlgebra.mul!(y::AbstractVector{T}, A::TimeDependentOperatorSum, x::AbstractVector, α, β) where T
# # Note that β is applied only to the first term
# mul!(y, A.operators[1], x, α*A.coefficients[1], β)
# @inbounds for i in 2:length(A.operators)
# mul!(y, A.operators[i], x, α*A.coefficients[i], 1)
# end
# y
# end

# @inline LinearAlgebra.mul!(y::AbstractVector{Ty}, A::QuantumObject{<:AbstractMatrix{Ta}}, x, α, β) where {Ty,Ta} = mul!(y, A.data, x, α, β)
## Sum of operators

mutable struct OperatorSum{CT<:Vector{<:Number},OT<:Vector{<:QuantumObject}}
coefficients::CT
Expand All @@ -107,6 +73,8 @@ end
Base.size(A::OperatorSum) = size(A.operators[1])
Base.size(A::OperatorSum, inds...) = size(A.operators[1], inds...)
Base.length(A::OperatorSum) = length(A.operators[1])
Base.copy(A::OperatorSum) = OperatorSum(copy(A.coefficients), copy(A.operators))
Base.deepcopy(A::OperatorSum) = OperatorSum(deepcopy(A.coefficients), deepcopy(A.operators))

function update_coefficients!(A::OperatorSum, coefficients)
length(A.coefficients) == length(coefficients) || throw(DimensionMismatch("The number of coefficients must be the same as the number of operators."))
Expand All @@ -117,6 +85,7 @@ end
# Note that β is applied only to the first term
mul!(y, A.operators[1], x, α*A.coefficients[1], β)
@inbounds for i in 2:length(A.operators)
A.coefficients[i] == 0 && continue
mul!(y, A.operators[i], x, α*A.coefficients[i], 1)
end
y
Expand Down
77 changes: 52 additions & 25 deletions src/time_evolution/time_evolution_dynamical.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,14 +232,17 @@ function _DSF_mesolve_Affect!(integrator)
dsf_params = internal_params.dsf_params
expv_cache = internal_params.expv_cache
dsf_identity = internal_params.dsf_identity
dsf_kron_cache_left = internal_params.dsf_kron_cache_left
dsf_kron_cache_left_dag = internal_params.dsf_kron_cache_left_dag
dsf_kron_cache_right = internal_params.dsf_kron_cache_right
dsf_kron_cache_right_dag = internal_params.dsf_kron_cache_right_dag
dsf_kron_cache_full = internal_params.dsf_kron_cache_full
# dsf_displace_cache_left = internal_params.dsf_displace_cache_left
# dsf_displace_cache_left_dag = internal_params.dsf_displace_cache_left_dag
# dsf_displace_cache_right = internal_params.dsf_displace_cache_right
# dsf_displace_cache_right_dag = internal_params.dsf_displace_cache_right_dag
dsf_displace_cache_full = internal_params.dsf_displace_cache_full

op_l_length = length(op_l)
dsf_displace_cache_full.coefficients .= 0

for i in eachindex(op_l)
op = op_l[i]
# op = op_l[i]
op_vec = op_l_vec[i]
αt = αt_list[i]
δα = δα_list[i]
Expand All @@ -252,17 +255,27 @@ function _DSF_mesolve_Affect!(integrator)

# This is equivalent to the code above, assuming that transpose(op) = adjoint(op)
# Aᵢ = kron(Δα * op.data - conj(Δα) * op.data', dsf_identity) + kron(dsf_identity, conj(Δα) * op.data - Δα * op.data')
@. dsf_kron_cache_full[i] = Δα * dsf_kron_cache_left[i] - conj(Δα) * dsf_kron_cache_left_dag[i] + conj(Δα) * dsf_kron_cache_right[i] - Δα * dsf_kron_cache_right_dag[i]
Aᵢ = dsf_kron_cache_full[i]

dsf_cache .= integrator.u
arnoldi!(expv_cache, Aᵢ, dsf_cache)
expv!(integrator.u, expv_cache, one(αt), dsf_cache)
# @. dsf_displace_cache_full[i] = Δα * dsf_displace_cache_left[i] - conj(Δα) * dsf_displace_cache_left_dag[i] + conj(Δα) * dsf_displace_cache_right[i] - Δα * dsf_displace_cache_right_dag[i]
# Aᵢ = dsf_displace_cache_full[i]

# dsf_cache .= integrator.u
# arnoldi!(expv_cache, Aᵢ, dsf_cache)
# expv!(integrator.u, expv_cache, one(αt), dsf_cache)

dsf_displace_cache_full.coefficients[i] = Δα
dsf_displace_cache_full.coefficients[i + op_l_length] = -conj(Δα)
dsf_displace_cache_full.coefficients[i + 2*op_l_length] = conj(Δα)
dsf_displace_cache_full.coefficients[i + 3*op_l_length] = -Δα

αt_list[i] += Δα
end
end

dsf_cache .= integrator.u
arnoldi!(expv_cache, dsf_displace_cache_full, dsf_cache)
expv!(integrator.u, expv_cache, 1, dsf_cache)

op_l2 = op_l .+ αt_list
e_ops2 = e_ops(op_l2, dsf_params)
_mat2vec_data = op -> mat2vec(get_data(op)')
Expand Down Expand Up @@ -294,20 +307,18 @@ function dsf_mesolveProblem(H::Function,
# Create the Krylov subspace with kron(H₀.data, H₀.data) just for initialize
expv_cache = arnoldi(kron(H₀.data, H₀.data), mat2vec(ket2dm(ψ0).data), krylov_dim)
dsf_identity = I(prod(H₀.dims))
dsf_kron_cache_left = map(op->kron(op.data, dsf_identity), op_l)
dsf_kron_cache_left_dag = map(op->kron(sparse(op.data'), dsf_identity), op_l)
dsf_kron_cache_right = map(op->kron(dsf_identity, op.data), op_l)
dsf_kron_cache_right_dag = map(op->kron(dsf_identity, sparse(op.data')), op_l)
dsf_kron_cache_full = @. dsf_kron_cache_left + dsf_kron_cache_left_dag + dsf_kron_cache_right + dsf_kron_cache_right_dag
dsf_displace_cache_left = map(op->Qobj(kron(op.data, dsf_identity)), op_l)
dsf_displace_cache_left_dag = map(op->Qobj(kron(sparse(op.data'), dsf_identity)), op_l)
dsf_displace_cache_right = map(op->Qobj(kron(dsf_identity, op.data)), op_l)
dsf_displace_cache_right_dag = map(op->Qobj(kron(dsf_identity, sparse(op.data'))), op_l)
dsf_displace_cache_full = OperatorSum(zeros(length(op_l) * 4), vcat(dsf_displace_cache_left, dsf_displace_cache_left_dag, dsf_displace_cache_right, dsf_displace_cache_right_dag))

params2 = params
params2 = merge(params, (H_fun = H, c_ops_fun = c_ops, e_ops_fun = e_ops,
op_l = op_l, op_l_vec = op_l_vec, αt_list = αt_list, δα_list = δα_list,
dsf_cache = similar(ψ0.data, length(ψ0.data)^2), expv_cache=expv_cache,
dsf_identity=dsf_identity, dsf_params = dsf_params,
dsf_kron_cache_left=dsf_kron_cache_left, dsf_kron_cache_left_dag=dsf_kron_cache_left_dag,
dsf_kron_cache_right=dsf_kron_cache_right, dsf_kron_cache_right_dag=dsf_kron_cache_right_dag,
dsf_kron_cache_full=dsf_kron_cache_full))
dsf_displace_cache_full=dsf_displace_cache_full))

cb_dsf = DiscreteCallback(_DSF_mesolve_Condition, _DSF_mesolve_Affect!, save_positions=(false, false))
kwargs2 = (;kwargs...)
Expand Down Expand Up @@ -412,6 +423,10 @@ function _DSF_mcsolve_Affect!(integrator)
dsf_cache = internal_params.dsf_cache2
expv_cache = internal_params.expv_cache
dsf_params = internal_params.dsf_params
dsf_displace_cache_full = internal_params.dsf_displace_cache_full

op_l_length = length(op_l)
dsf_displace_cache_full.coefficients .= 0

for i in eachindex(op_l)
op = op_l[i]
Expand All @@ -424,16 +439,22 @@ function _DSF_mcsolve_Affect!(integrator)
# dsf_cache .= integrator.u
# mul!(integrator.u, Dᵢ.data', dsf_cache)

Aᵢ = -Δα*op.data' + conj(Δα)*op.data
dsf_cache .= integrator.u
arnoldi!(expv_cache, Aᵢ, dsf_cache)
expv!(integrator.u, expv_cache, one(αt), dsf_cache)
# Aᵢ = -Δα*op.data' + conj(Δα)*op.data
# dsf_cache .= integrator.u
# arnoldi!(expv_cache, Aᵢ, dsf_cache)
# expv!(integrator.u, expv_cache, one(αt), dsf_cache)

dsf_displace_cache_full.coefficients[i] = conj(Δα)
dsf_displace_cache_full.coefficients[i + op_l_length] = -Δα

αt_list[i] += Δα
end
end

dsf_cache .= integrator.u
arnoldi!(expv_cache, dsf_displace_cache_full, dsf_cache)
expv!(integrator.u, expv_cache, 1, dsf_cache)

op_l2 = op_l .+ αt_list
e_ops2 = e_ops(op_l2, dsf_params)
c_ops2 = c_ops(op_l2, dsf_params)
Expand All @@ -452,7 +473,9 @@ function _dsf_mcsolve_prob_func(prob, i, repeat)
cumsum_weights_mc = similar(internal_params.weights_mc), random_n = Ref(rand()), progr_mc = ODEProgress(0),
jump_times = similar(internal_params.jump_times), jump_which = similar(internal_params.jump_which),
αt_list = deepcopy(internal_params.αt_list), dsf_cache1 = similar(internal_params.dsf_cache1),
dsf_cache2 = similar(internal_params.dsf_cache2), expv_cache = deepcopy(internal_params.expv_cache)))
dsf_cache2 = similar(internal_params.dsf_cache2), expv_cache = deepcopy(internal_params.expv_cache),
# dsf_displace_cache_full = deepcopy(internal_params.dsf_displace_cache_full),
dsf_displace_cache_full = OperatorSum(deepcopy(internal_params.dsf_displace_cache_full.coefficients), internal_params.dsf_displace_cache_full.operators)))

remake(prob, p=prm)
end
Expand Down Expand Up @@ -480,11 +503,15 @@ function dsf_mcsolveEnsembleProblem(H::Function,
αt_list = convert(Vector{T}, α0_l)
expv_cache = arnoldi(H₀.data, ψ0.data, krylov_dim)

dsf_displace_cache = map(op->Qobj(op.data), op_l)
dsf_displace_cache_dag = map(op->Qobj(sparse(op.data')), op_l)
dsf_displace_cache_full = OperatorSum(zeros(length(op_l) * 2), vcat(dsf_displace_cache, dsf_displace_cache_dag))


params2 = merge(params, (H_fun = H, c_ops_fun = c_ops, e_ops_fun = e_ops,
op_l = op_l, αt_list = αt_list, δα_list = δα_list,
dsf_cache1 = similar(ψ0.data), dsf_cache2 = similar(ψ0.data),
expv_cache = expv_cache, dsf_params=dsf_params))
expv_cache = expv_cache, dsf_params=dsf_params, dsf_displace_cache_full=dsf_displace_cache_full))

cb_dsf = DiscreteCallback(_DSF_mcsolve_Condition, _DSF_mcsolve_Affect!, save_positions=(false, false))
kwargs2 = (;kwargs...)
Expand Down

0 comments on commit 3c53b1a

Please sign in to comment.