|
1 | | -export mesolveProblem, mesolve |
| 1 | +export mesolveProblem, mesolve, mesolve_map |
2 | 2 |
|
3 | 3 | _mesolve_make_L_QobjEvo(H::Union{QuantumObject,Nothing}, c_ops) = QobjEvo(liouvillian(H, c_ops); type = SuperOperator()) |
4 | 4 | _mesolve_make_L_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}, c_ops) = liouvillian(QobjEvo(H), c_ops) |
5 | 5 | _mesolve_make_L_QobjEvo(H::Nothing, c_ops::Nothing) = throw(ArgumentError("Both H and |
6 | 6 | c_ops are Nothing. You are probably running the wrong function.")) |
7 | 7 |
|
| 8 | +function _gen_mesolve_solution(sol, times, dimensions, isoperket::Val) |
| 9 | + if getVal(isoperket) |
| 10 | + ρt = map(ϕ -> QuantumObject(ϕ, type = OperatorKet(), dims = dimensions), sol.u) |
| 11 | + else |
| 12 | + ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator(), dims = dimensions), sol.u) |
| 13 | + end |
| 14 | + |
| 15 | + kwargs = NamedTuple(sol.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility |
| 16 | + |
| 17 | + return TimeEvolutionSol( |
| 18 | + times, |
| 19 | + sol.t, |
| 20 | + ρt, |
| 21 | + _get_expvals(sol, SaveFuncMESolve), |
| 22 | + sol.retcode, |
| 23 | + sol.alg, |
| 24 | + kwargs.abstol, |
| 25 | + kwargs.reltol, |
| 26 | + ) |
| 27 | +end |
| 28 | + |
8 | 29 | @doc raw""" |
9 | 30 | mesolveProblem( |
10 | 31 | H::Union{AbstractQuantumObject,Tuple}, |
@@ -207,23 +228,143 @@ end |
207 | 228 | function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); kwargs...) |
208 | 229 | sol = solve(prob.prob, alg; kwargs...) |
209 | 230 |
|
210 | | - # No type instabilities since `isoperket` is a Val, and so it is known at compile time |
211 | | - if getVal(prob.kwargs.isoperket) |
212 | | - ρt = map(ϕ -> QuantumObject(ϕ, type = OperatorKet(), dims = prob.dimensions), sol.u) |
213 | | - else |
214 | | - ρt = map(ϕ -> QuantumObject(vec2mat(ϕ), type = Operator(), dims = prob.dimensions), sol.u) |
| 231 | + return _gen_mesolve_solution(sol, prob.times, prob.dimensions, prob.kwargs.isoperket) |
| 232 | +end |
| 233 | + |
| 234 | +@doc raw""" |
| 235 | + mesolve_map( |
| 236 | + H::Union{AbstractQuantumObject,Tuple}, |
| 237 | + ψ0::Union{QuantumObject,AbstractVector{<:QuantumObject}}, |
| 238 | + tlist::AbstractVector, |
| 239 | + c_ops::Union{Nothing,AbstractVector,Tuple} = nothing; |
| 240 | + alg::OrdinaryDiffEqAlgorithm = Tsit5(), |
| 241 | + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), |
| 242 | + e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, |
| 243 | + params::Union{NullParameters,Tuple} = NullParameters(), |
| 244 | + progress_bar::Union{Val,Bool} = Val(true), |
| 245 | + kwargs..., |
| 246 | + ) |
| 247 | +
|
| 248 | +Solve the master equation for multiple initial states and parameter sets using ensemble simulation. |
| 249 | +
|
| 250 | +This function computes the time evolution for all combinations (Cartesian product) of initial states and parameter sets, solving the Lindblad master equation (see [`mesolve`](@ref)): |
| 251 | +
|
| 252 | +```math |
| 253 | +\frac{\partial \hat{\rho}(t)}{\partial t} = -i[\hat{H}, \hat{\rho}(t)] + \sum_n \mathcal{D}(\hat{C}_n) [\hat{\rho}(t)] |
| 254 | +``` |
| 255 | +
|
| 256 | +where |
| 257 | +
|
| 258 | +```math |
| 259 | +\mathcal{D}(\hat{C}_n) [\hat{\rho}(t)] = \hat{C}_n \hat{\rho}(t) \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n \hat{\rho}(t) - \frac{1}{2} \hat{\rho}(t) \hat{C}_n^\dagger \hat{C}_n |
| 260 | +``` |
| 261 | +
|
| 262 | +for each combination in the ensemble. |
| 263 | +
|
| 264 | +# Arguments |
| 265 | +
|
| 266 | +- `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs. |
| 267 | +- `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref). |
| 268 | +- `tlist`: List of time points at which to save either the state or the expectation values of the system. |
| 269 | +- `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`. |
| 270 | +- `alg`: The algorithm for the ODE solver. The default is `Tsit5()`. |
| 271 | +- `ensemblealg`: Ensemble algorithm to use for parallel computation. Default is `EnsembleThreads()`. |
| 272 | +- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`. |
| 273 | +- `params`: A `Tuple` of parameter sets. Each element should be an `AbstractVector` representing the sweep range for that parameter. The function will solve for all combinations of initial states and parameter sets. |
| 274 | +- `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities. |
| 275 | +- `kwargs`: The keyword arguments for the ODEProblem. |
| 276 | +
|
| 277 | +# Notes |
| 278 | +
|
| 279 | +- The function returns an array of solutions with dimensions matching the Cartesian product of initial states and parameter sets. |
| 280 | +- If `ψ0` is a vector of `m` states and `params = (p1, p2, ...)` where `p1` has length `n1`, `p2` has length `n2`, etc., the output will be of size `(m, n1, n2, ...)`. |
| 281 | +- If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolve_map`](@ref) instead. |
| 282 | +- See [`mesolve`](@ref) for more details. |
| 283 | +
|
| 284 | +# Returns |
| 285 | +
|
| 286 | +- An array of [`TimeEvolutionSol`](@ref) objects with dimensions `(length(ψ0), length(params[1]), length(params[2]), ...)`. |
| 287 | +""" |
| 288 | +function mesolve_map( |
| 289 | + H::Union{AbstractQuantumObject{HOpType},Tuple}, |
| 290 | + ψ0::AbstractVector{<:QuantumObject{StateOpType}}, |
| 291 | + tlist::AbstractVector, |
| 292 | + c_ops::Union{Nothing,AbstractVector,Tuple} = nothing; |
| 293 | + alg::OrdinaryDiffEqAlgorithm = Tsit5(), |
| 294 | + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), |
| 295 | + e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, |
| 296 | + params::Union{NullParameters,Tuple} = NullParameters(), |
| 297 | + progress_bar::Union{Val,Bool} = Val(true), |
| 298 | + kwargs..., |
| 299 | +) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}} |
| 300 | + (isoper(H) && all(isket, ψ0) && isnothing(c_ops)) && return sesolve_map( |
| 301 | + H, |
| 302 | + ψ0, |
| 303 | + tlist; |
| 304 | + alg = alg, |
| 305 | + ensemblealg = ensemblealg, |
| 306 | + e_ops = e_ops, |
| 307 | + params = params, |
| 308 | + progress_bar = progress_bar, |
| 309 | + kwargs..., |
| 310 | + ) |
| 311 | + |
| 312 | + # mapping initial states and parameters |
| 313 | + # Convert to appropriate format based on state type |
| 314 | + ψ0_iter = map(ψ0) do state |
| 315 | + T = _complex_float_type(eltype(state)) |
| 316 | + if isoperket(state) |
| 317 | + to_dense(T, copy(state.data)) |
| 318 | + else |
| 319 | + to_dense(T, mat2vec(ket2dm(state).data)) |
| 320 | + end |
215 | 321 | end |
| 322 | + iter = |
| 323 | + params isa NullParameters ? collect(Iterators.product(ψ0_iter, [params])) : |
| 324 | + collect(Iterators.product(ψ0_iter, params...)) |
| 325 | + ntraj = length(iter) |
216 | 326 |
|
217 | | - kwargs = NamedTuple(sol.prob.kwargs) # Convert to NamedTuple for Zygote.jl compatibility |
| 327 | + # we disable the progress bar of the mesolveProblem because we use a global progress bar for all the trajectories |
| 328 | + prob = mesolveProblem( |
| 329 | + H, |
| 330 | + first(ψ0), |
| 331 | + tlist, |
| 332 | + c_ops; |
| 333 | + e_ops = e_ops, |
| 334 | + params = first(iter)[2:end], |
| 335 | + progress_bar = Val(false), |
| 336 | + kwargs..., |
| 337 | + ) |
218 | 338 |
|
219 | | - return TimeEvolutionSol( |
| 339 | + # generate and solve ensemble problem |
| 340 | + _output_func = _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _standard_output_func) # setup global progress bar |
| 341 | + ens_prob = TimeEvolutionProblem( |
| 342 | + EnsembleProblem( |
| 343 | + prob.prob, |
| 344 | + prob_func = (prob, i, repeat) -> _se_me_map_prob_func(prob, i, repeat, iter), |
| 345 | + output_func = _output_func[1], |
| 346 | + safetycopy = false, |
| 347 | + ), |
220 | 348 | prob.times, |
221 | | - sol.t, |
222 | | - ρt, |
223 | | - _get_expvals(sol, SaveFuncMESolve), |
224 | | - sol.retcode, |
225 | | - sol.alg, |
226 | | - kwargs.abstol, |
227 | | - kwargs.reltol, |
| 349 | + prob.dimensions, |
| 350 | + (progr = _output_func[2], channel = _output_func[3], isoperket = prob.kwargs.isoperket), |
228 | 351 | ) |
| 352 | + sol = _ensemble_dispatch_solve(ens_prob, alg, ensemblealg, ntraj) |
| 353 | + |
| 354 | + # handle solution and make it become an Array of TimeEvolutionSol |
| 355 | + sol_vec = |
| 356 | + [_gen_mesolve_solution(sol[:, i], prob.times, prob.dimensions, prob.kwargs.isoperket) for i in eachindex(sol)] # map is type unstable |
| 357 | + if params isa NullParameters # if no parameters specified, just return a Vector |
| 358 | + return sol_vec |
| 359 | + else |
| 360 | + return reshape(sol_vec, size(iter)) |
| 361 | + end |
229 | 362 | end |
| 363 | +mesolve_map( |
| 364 | + H::Union{AbstractQuantumObject{HOpType},Tuple}, |
| 365 | + ψ0::QuantumObject{StateOpType}, |
| 366 | + tlist::AbstractVector, |
| 367 | + c_ops::Union{Nothing,AbstractVector,Tuple} = nothing; |
| 368 | + kwargs..., |
| 369 | +) where {HOpType<:Union{Operator,SuperOperator},StateOpType<:Union{Ket,Operator,OperatorKet}} = |
| 370 | + mesolve_map(H, [ψ0], tlist, c_ops; kwargs...) |
0 commit comments