Skip to content

Commit ee5f637

Browse files
committed
Implement support for OrdinaryDiffEq
This also includes a proper definition of "parametrized controls" with a `QuantumPropgagators.Controls.get_parameters` routine and a `QuantumPropgagators.Controls.ParameterizedFunction` abstract type.
1 parent e62fce4 commit ee5f637

25 files changed

+987
-90
lines changed

Project.toml

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,14 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1515
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
1616
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
1717

18+
[weakdeps]
19+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
20+
21+
[extensions]
22+
QuantumPropagatorsODEExt = "OrdinaryDiffEq"
23+
1824
[compat]
25+
OrdinaryDiffEq = "6.59"
1926
OffsetArrays = "1"
2027
ProgressMeter = "1"
2128
SpecialFunctions = "1, 2"

devrepl.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,10 @@ if !isfile(joinpath("test", "Manifest.toml"))
4242
end
4343
include("test/init.jl")
4444

45+
# Disable link-checking in interactive REPL, since it is the slowest part
46+
# of building the docs.
47+
ENV["DOCUMENTER_CHECK_LINKS"] = "0"
48+
4549
if abspath(PROGRAM_FILE) == @__FILE__
4650
help()
4751
end

docs/make.jl

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
using QuantumPropagators
2+
using QuantumPropagators: AbstractPropagator, set_t!, set_state!
23
using Documenter
34
using Pkg
5+
using OrdinaryDiffEq # ensure ODE extension is loaded
46
using DocumenterCitations
57

68
DocMeta.setdocmeta!(
@@ -22,10 +24,19 @@ include("generate_api.jl")
2224

2325
bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib"); style=:numeric)
2426

27+
warnonly = [:linkcheck,]
28+
if get(ENV, "DOCUMENTER_WARN_ONLY", "0") == "1" # cf. test/init.jl
29+
warnonly = true
30+
end
31+
32+
2533
makedocs(;
2634
plugins=[bib],
2735
authors=AUTHORS,
2836
sitename="QuantumPropagators.jl",
37+
# Link checking is disabled in REPL, see `devrepl.jl`.
38+
linkcheck=(get(ENV, "DOCUMENTER_CHECK_LINKS", "1") != "0"),
39+
warnonly,
2940
format=Documenter.HTML(;
3041
prettyurls=true,
3142
canonical="https://juliaquantumcontrol.github.io/QuantumPropagators.jl",
@@ -38,10 +49,6 @@ makedocs(;
3849
"Dynamical Generators" => "generators.md",
3950
"Propagation Methods" => "methods.md",
4051
"Expectation Values" => "storage.md",
41-
hide("Examples" => joinpath("examples", "index.md"), [
42-
# "Example 1" => joinpath("examples", "1.md"),
43-
# "Example 2" => joinpath("examples", "2.md"),
44-
]),
4552
"Howtos" => "howto.md",
4653
hide("Benchmarks" => "benchmarks.md", [joinpath("benchmarks", "profiling.md")]),
4754
"API" => "api/quantumpropagators.md",

docs/src/benchmarks/profiling.md

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,3 @@
1-
```@meta
2-
CurrentModule = QuantumPropagators
3-
```
4-
51
# Profiling Howto
62

73
To choose an appropriate propagation method and parameters for a given problem, it is essential to benchmark and profile the propagation.

docs/src/examples/index.md

Lines changed: 0 additions & 12 deletions
This file was deleted.

docs/src/howto.md

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,19 @@
11
# Howtos
22

3-
## How to extend QuantumPropagators with a new propagation method
3+
## How to implement a new propagation method
44

55
* Define a new sub-type of [`AbstractPropagator`](@ref QuantumPropagators.AbstractPropagator) type that is unique to the propagation method, e.g. `MyNewMethodPropagator`. If appropriate, sub-type [`PiecewisePropagator`](@ref QuantumPropagators.PiecewisePropagator) or [`PWCPropagator`](@ref QuantumPropagators.PWCPropagator).
66

7-
* Choose a name for the propagation method, e.g. `mynewmethod` and implement a new [`init_prop`](@ref) method with the following signature
7+
* The high-level [`propagate`](@ref) and [`init_prop`](@ref) functions have a mandatory `method` keyword argument. That argument should receive a `Module` object for the `module` implementing a propagation method (e.g., `using QuantumPropagators: Cheby; method=Cheby`). This ensures that the module or package implementing the method is fully loaded. Internally, [`init_prop`](@ref) delegates `method=module` to a *positional* argument `method = Val(nameof(module))`
8+
9+
* Thus, if `MyNewMethodPropagator` is implemented in a module `MyNewMethod`, or if it wraps a registered package `MyNewMethod`, implement a new [`init_prop`](@ref) method with the following signature:
810

911
```
1012
function init_prop(
1113
state,
1214
generator,
1315
tlist,
14-
method::Val{:mynewmethod};
16+
method::Val{:MyNewMethod};
1517
inplace=true,
1618
backward=false,
1719
verbose=false,
@@ -21,6 +23,10 @@
2123
)
2224
```
2325

24-
Note the `method::Val{:mynewmethod}` as the fourth positional parameter. While the *public* interface for [`init_prop`](@ref) takes `method` as a keyword argument, privately [`init_prop`](@ref) dispatches for different methods as above.
26+
Note the `method::Val{:MyNewMethod}` as the fourth positional (!) parameter. While the *public* interface for [`init_prop`](@ref) takes `method` as a keyword argument, privately [`init_prop`](@ref) dispatches for different methods as above.
27+
28+
* If the propagation method is not associated with a module or package, it is also possible to implement a method [`init_prop`](@ref) with a fourth positional argument with a type of, e.g., `::Val{:my_new_method}`. This would allow to call the high-level [`propagate`](@ref)/[`init_prop`](@ref) with the keyword argument `method=:my_new_method`, i.e., passing a name (`Symbol`) instead of a `Module`.
29+
30+
* Implement the remaining methods in [The Propagator interface](@ref).
2531

26-
* Implement the remaining methods in [The Propagator interface](@ref)
32+
* Test the implementation by instantiating a `propagator` and calling [`QuantumPropagators.Interfaces.check_propagator`](@ref) on it.

docs/src/index.md

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,3 @@
1-
```@meta
2-
CurrentModule = QuantumPropagators
3-
```
4-
51
# QuantumPropagators
62

73
```@eval
@@ -19,6 +15,14 @@ Markdown.parse("$github_badge $version_badge")
1915

2016
The [QuantumPropagators](https://github.com/JuliaQuantumControl/QuantumPropagators.jl) package implements methods for simulating the time dynamics of a quantum system. It is the numerical backend for all packages implementing methods of optimal control within the [JuliaQuantumControl](https://github.com/JuliaQuantumControl) organization.
2117

18+
The purpose is this package is to:
19+
20+
* Implement the piecewise-constant propagation methods that have been traditionally used in quantum control, based on an expansion of the time evolution operator into Chebychev polynomials (closed systems) or Newton polynomials (open systems),
21+
* Delegate to third-party packages that implement state-of-the-art methods for time-continuous dynamics ([DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/)).
22+
* Define data structures for "dynamical generators" (time-dependent Hamiltonians or Liouvillians) that are relevant to quantum control. Specifically, this includes a well-defined notion of "control functions" and "control parameters" inside the Hamiltonian.
23+
* Define a general interface for a time propagators that is suitable for all aspects of quantum control. This includes the ability to efficiently tune control parameters or control functions while the propagation is running (for "sequential updates"), and the ability to propagate both forwards and backwards in time.
24+
25+
The focus on quantum control separates this package from other solutions for just simulating the time dynamics of a quantum system, e.g. within [QuantumOptics.jl](https://qojulia.org) or directly with [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/). Also, while `QuantumPropagators` defines a general interface for dynamic generators, it is otherwise agnostic about the data structures to represent Hamiltonians or quantum states. This allows for the use of modeling frameworks like [QuantumOptics.jl](https://qojulia.org) or of custom problem-specific data structures.
2226

2327

2428
## Contents
@@ -29,14 +33,14 @@ Pages=[
2933
"generators.md",
3034
"methods.md",
3135
"storage.md",
32-
"examples/index.md",
3336
"howto.md",
3437
"benchmarks.md",
3538
"api/quantumpropagators.md",
3639
]
3740
Depth = 2
3841
```
3942

43+
4044
## History
4145

4246
See the [Releases](https://github.com/JuliaQuantumControl/QuantumPropagators.jl/releases) on Github.

docs/src/methods.md

Lines changed: 159 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,168 @@
11
# Propagation Methods
22

3-
## Chebychev Propagation
3+
```math
4+
\gdef\op#1{\hat{#1}}
5+
\gdef\ket#1{\vert{#1}\rangle}
6+
\gdef\Liouvillian{\mathcal{L}}
7+
\gdef\Re{\operatorname{Re}}
8+
\gdef\Im{\operatorname{Im}}
9+
```
410

5-
For `method=:cheby`, the time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Chebychev polynomials [Tal-EzerJCP1984](@cite)[KosloffJCP1988](@cite). This requires ```` to be Hermitian (have real eigenvalues) and to have a known spectral range, so that it can be normalized to the domain ``[-1, 1]`` on which the Chebychev polynomials are defined.
11+
As discussed in the [Overview](@ref overview_approaches), time propagation can be implemented in one of two ways:
612

7-
… (TODO) …
13+
1. By *analytically* solving the equation of motion and numerically evaluating the application time evolution operator.
814

9-
## Newton Propagation
15+
We consider this especially in the piecewise-constant case (`pwc=true` in [`propagate`](@ref)/[`init_prop`](@ref)), which is required for the traditional optimization methods [GRAPE](https://juliaquantumcontrol.github.io/GRAPE.jl/stable/) and [Krotov](https://juliaquantumcontrol.github.io/Krotov.jl/stable/). In these propagations, the time-dependent generator ``\op{H}(t)`` is [evaluated](@ref QuantumPropagators.Controls.evaluate) to a constant operator ``\op{H}`` on each interval of the time grid. The analytical solution to the Schrödinger or Liouville equation is well known, and propagation step simply has to evaluate the application of the time evolution operator ``\op{U} = \exp[-i \op{H} dt]`` to the state ``|Ψ⟩``. The following methods are built in to `QuantumPropagators`:
1016

11-
For `method=:newton`, the time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Newton polynomials [BermanJPA1992](@cite)[AshkenaziJCP1995](@cite)[Tal-EzerSJSC2007](@cite). Unlike for Chebychev polynomials, this expansion does not require ```` to be Hermitian or to have a known spectral radius. This makes the Newton propagation applicable to open quantum systems, where ```` is replaced by a Liouvillian to calculate the time evolution of the density matrix.
17+
* [`ExpProp`](@ref method_expprop) – constructs ``\op{U}`` explicitly and then applies it to ``|Ψ⟩``
18+
* [`Cheby`](@ref method_cheby) — expansion of ``\op{U} |Ψ⟩`` into Chebychev polynomials, valid if ``\op{H}`` has real eigenvalues
19+
* [`Newton`](@ref method_newton) – expansion of ``\op{U} |Ψ⟩`` into Newton polynomials, valid if ``\op{H}`` has complex eigenvalues (non-Hermitian Hamiltonian, Liouvillian)
1220

13-
… (TODO) …
21+
The `ExpProp` method is generally not numerically efficient, but works well for small system for for debugging. The two "core" methods based on a polynomials series expansions are more suitable for bigger systems and provide both efficiency and high precision (in general, the is truncated as soon as some desired precision is reached, which is machine precision by default).
1422

15-
## Propagation with Explicit Matrix Exponentiation
23+
Note that this high precision is *within the piecewise-constant approximation*. The discretization itself may introduce a non-negligible error compared to the time-continuous dynamics. There is tradeoff: A smaller `dt` decreases the discretization error, but the polynomial expansions are more effective with larger time steps.
1624

17-
TODO
25+
2. By solving the equation of motion explicitly with an ODE solver.
26+
27+
We support the use of any of the ODE solvers in [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/):
28+
29+
* [`OrdinaryDiffEq`](@ref method_ode) – solve the equation of motion as an ODE
30+
31+
The main benefit of using an ODE solver is that the generator ``\op{H}(t)`` can be treated as time-continuous, and thus avoid the time discretization error. While this is not compatible with traditional optimal control method like [GRAPE](https://juliaquantumcontrol.github.io/GRAPE.jl/stable/) and [Krotov](https://juliaquantumcontrol.github.io/Krotov.jl/stable/), it is suitable for control methods for tuning analytical control parameters [LucarelliPRA2018,MachnesPRL2018,SorensenPRA2018](@cite).
32+
33+
The `method=OrdinaryDiffEq` is also available in a piecewise-constant mode by setting `pwc=true`, for comparison with `method=Cheby` and `method=Newton`.
34+
35+
36+
## [`ExpProp`](@id method_expprop)
37+
38+
The method should be loaded with
39+
40+
```
41+
using QuantumPropagators: ExpProp
42+
```
43+
44+
and the passed as `method=ExpProp` to [`propagate`](@ref) or [`init_prop`](@ref):
45+
46+
```@docs
47+
init_prop(state, generator, tlist, method::Val{:ExpProp}; kwargs...)
48+
```
49+
50+
**Advantages**
51+
52+
* Simple: no knobs to turn
53+
* "Exact" to machine precision (within the piecewise constant approximation)
54+
* Does not require any special properties or knowledge of the dynamical generator
55+
* Efficient for small systems
56+
57+
**Disadvantages**
58+
59+
* Bad numerical scaling with the Hilbert space dimension
60+
* Method for `exp(-1im * H * dt)` must be defined (or `H` must be convertible to a type that can be exponentiated)
61+
62+
**When to use**
63+
64+
* Small Hilbert space dimension (<10)
65+
* Comparing against another propagator
66+
67+
68+
## [`Cheby`](@id method_cheby)
69+
70+
The method should be loaded with
71+
72+
```
73+
using QuantumPropagators: Cheby
74+
```
75+
76+
and then passed as `method=Cheby` to [`propagate`](@ref) or [`init_prop`](@ref):
77+
78+
```@docs
79+
init_prop(state, generator, tlist, method::Val{:Cheby}; kwargs...)
80+
```
81+
82+
The time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Chebychev polynomials [Tal-EzerJCP1984, KosloffJCP1988](@cite). This requires ```` to be Hermitian (have real eigenvalues) and to have a known spectral range, so that it can be normalized to the domain ``[-1, 1]`` on which the Chebychev polynomials are defined.
83+
84+
See [GoerzPhd2015; Chapter 3.2.1](@cite) for a detailed description of the method.
85+
86+
**Advantages**
87+
88+
* Very efficient for high precision and moderately large time steps
89+
90+
**Disadvantages**
91+
92+
* Only valid for Hermitian operators
93+
* Must be able to estimate the spectral envelope
94+
95+
**When to use**
96+
97+
* Closed quantum systems with piecewise constant dynamics
98+
99+
100+
## [`Newton`](@id method_newton)
101+
102+
The method should be loaded with
103+
104+
```
105+
using QuantumPropagators: Newton
106+
```
107+
108+
and then passed as `method=Newton` to [`propagate`](@ref) or [`init_prop`](@ref):
109+
110+
111+
```@docs
112+
init_prop(state, generator, tlist, method::Val{:Newton}; kwargs...)
113+
```
114+
115+
The time evolution operator of the piecewise-constant Schrödinger equation ``|Ψ(t)⟩ = e^{-i Ĥ dt} |Ψ(0)⟩`` is evaluated by an expansion into Newton polynomials [BermanJPA1992, AshkenaziJCP1995, Tal-EzerSJSC2007](@cite). Unlike for Chebychev polynomials, this expansion does not require ```` to be Hermitian or to have a known spectral radius. This makes the Newton propagation applicable to open quantum systems, where ```` is replaced by a Liouvillian to calculate the time evolution of the density matrix.
116+
117+
See [GoerzPhd2015; Chapter 3.2.2](@cite) for a detailed description of the method.
118+
119+
**Advantages**
120+
121+
* Reasonably efficient for high precision and moderately large time steps
122+
* Spectral radius does not need to be known
123+
124+
**Disadvantages**
125+
126+
* Need to choose `m_max` and `max_restarts` well for good performance.
127+
128+
**When to use**
129+
130+
* Open quantum systems with piecewise constant dynamics
131+
132+
133+
## [`OrdinaryDiffEq`](@id method_ode)
134+
135+
The method requires that the [OrdinaryDiffEq](https://docs.sciml.ai/OrdinaryDiffEq/stable/) package is loaded
136+
137+
```
138+
using OrdinaryDiffEq
139+
```
140+
141+
Equivalently, the more general [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package can be used.
142+
143+
```
144+
using DifferentialEquations
145+
```
146+
147+
There is no difference between these two options: `OrdinaryDiffEq` is just a smaller dependency, but `DifferentialEquations` may be preferred if the large DifferentialEquations framework is required for the project.
148+
149+
In any case, the loaded package to [`propagate`](@ref) or [`init_prop`](@ref) via the `method` keyword argument:
150+
151+
152+
```@docs
153+
init_prop(state, generator, tlist, method::Val{:OrdinaryDiffEq}; kwargs...)
154+
```
155+
156+
**Advantages**
157+
158+
* Suitable for time-continuous dynamics
159+
* The full power of the DifferentialEquations.jl ecosystem
160+
* Efficient for moderate precisions
161+
162+
**Disadvantages**
163+
164+
* Less efficient for piecewise-constant dynamics, and thus less suitable of PWC control methods
165+
166+
**When to use**
167+
168+
* Time-continuous dynamics

0 commit comments

Comments
 (0)