CTFlows.jl

Appearance

Flow(HamiltonianVectorField) compatibility ​

This page is a living compatibility reference for the Hamiltonian flow built from a Data.HamiltonianVectorField: which state/costate types and call styles it accepts, each shown with a minimal, executable example. Every ✓ / ⚠ in the table below is demonstrated by a code block on this page and is re-run on every documentation build, so the page cannot drift from the code.

Scope: CPU only. The table is generated from probe/cpu/probe_cpu.jl (run it locally with julia --project=probe/cpu probe/cpu/probe_cpu.jl). GPU compatibility is a separate effort — see probe/gpu.

All examples integrate the harmonic oscillator     (autonomous, fixed) with the default SciML integrator (OrdinaryDiffEqTsit5). The analytic solution is   ,    , so integrating from   to   maps every initial condition to   .

Last probed: 2026-07-23

This page's table shape — which containers and axes are tested — reflects probe/cpu/probe_cpu.jl as of this date. Every ✓/⚠ cell below is still re-executed on every documentation build regardless (see Compatibility overview) — only the scope of what's tested can go stale, not the results shown.

Two flows are built once — one from an out-of-place Hamiltonian vector field, one from an in-place one — and reused by every example below, so each example isolates the one thing that varies: the initial condition (x0, p0).

julia
hvf   = Data.HamiltonianVectorField((x, p) -> (p, -x))   # out-of-place: (x,p) -> (p,-x)
hflow = Flows.Flow(hvf; reltol=1e-8)
Flow
├─ system: HamiltonianVectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  └─ HamiltonianVectorField: autonomous, fixed (no variable), out-of-place
│       natural call: f(x, p)
│       uniform call: f(t, x, p, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.
julia
hvf_ip   = Data.HamiltonianVectorField(
    (dx, dp, x, p) -> (dx .= p; dp .= -x); is_autonomous=true, is_variable=false,
)   # in-place
hflow_ip = Flows.Flow(hvf_ip; reltol=1e-8)
Flow
├─ system: HamiltonianVectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  └─ HamiltonianVectorField: autonomous, fixed (no variable), in-place
│       natural call: f(dx, dp, x, p)
│       uniform call: f(dx, dp, t, x, p, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.

Compatibility table ​

The two call styles are the point call hflow(t0, x0, p0, tf) → final (xf, pf), and the trajectory call hflow((t0, tf), x0, p0) → Trajectories.HamiltonianVectorFieldTrajectory. The OOP / IP columns are the Hamiltonian-vector-field kind — an out-of-place (x, p) -> (p, -x) versus an in-place (dx, dp, x, p) -> (dx .= p; dp .= -x) — independent of the state/costate container. (Which integration path is taken is chosen automatically from the mutability of vcat(x0, p0); see In-place Hamiltonian vector fields.)

State/costate typeOOP pointOOP trajIP pointIP traj
Scalar Real✓✓✓✓
Vector Real✓✓✓✓
MVector Real✓✓✓✓
SVector Real✓✓⚠⚠
Matrix Real (batch)✓✓✓✓
MMatrix Real✓✓✓✓
SMatrix Real✓✓⚠⚠
Scalar Complex✓✓✓✓
Vector Complex✓✓✓✓
MVector Complex✓✓✓✓
SVector Complex✓✓⚠⚠
Matrix Complex✓✓✓✓
SMatrix Complex✓✓⚠⚠
ForwardDiff.Dual scalar✓✓✓✓
ForwardDiff.Dual Vector✓✓✓✓
ForwardDiff.Dual MVector✓✓✓✓
ForwardDiff.Dual SVector✓✓⚠⚠

Legend ​

  • ✓ — works; the result matches the analytic    and is executed on this page.

  • ⚠ — works, but emits a performance warning; see note (a).

Every state/costate type is supported on CPU — there is no unsupported (✗) combination.

(a) In-place Hamiltonian vector field + immutable initial condition

An in-place Hamiltonian vector field with an immutable (x0, p0) (SVector, SMatrix) emits a performance @warn and falls back to an out-of-place finalize RHS. The result is correct but the path is slower. Prefer an out-of-place Hamiltonian vector field, or mutable MVector / MMatrix, for static states.

1-D = scalar, end to end

A scalar (x0, p0) stays a scalar in both call styles: the point call returns a scalar (xf, pf) pair, and the trajectory accessors state(sol)(t) / costate(sol)(t) also return scalars — not length-1 vectors. A length-1 vector or SVector collapses the same way. This is the ecosystem's "1-D = scalar" convention, implemented here via the HamiltonianVectorFieldSystem coercion machinery — the same mechanism Flow(VectorField) now shares. See the Shape contract page for the full cross-constructor picture.

Real states ​

Scalar ​

Both call styles return scalars:

julia
julia> hflow(0.0, 1.0, 0.0, pi/2)   # (xf, pf) ≈ (0.0, -1.0)
(-1.5541996200011742e-10, -0.9999999997009266)
julia
julia> sol = hflow((0.0, pi/2), 1.0, 0.0);

julia> (Trajectories.state(sol)(pi/2), Trajectories.costate(sol)(pi/2))   # ≈ (0.0, -1.0), both scalars
(-1.5541963451763603e-10, -0.9999999997009266)

Vector ​

julia
julia> hflow(0.0, [1.0, 0.0], [0.0, 1.0], pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

The trajectory's state/costate accessors interpolate at any time in the span:

julia
julia> sol = hflow((0.0, pi/2), [1.0, 0.0], [0.0, 1.0]);

julia> (Trajectories.state(sol)(pi/4), Trajectories.costate(sol)(pi/4))
([0.7071067807847535, 0.707106781788745], [-0.707106781788745, 0.7071067807847535])

MVector and SVector ​

Mutable and immutable static vectors both integrate out-of-place:

julia
julia> hflow(0.0, MVector{2}(1.0, 0.0), MVector{2}(0.0, 1.0), pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

julia> hflow(0.0, SA[1.0, 0.0], SA[0.0, 1.0], pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

Matrix (batch) ​

A matrix state/costate integrates every column as an independent trajectory of the harmonic oscillator:

julia
julia> hflow(0.0, [1.0 2.0; 3.0 4.0], [0.0 0.0; 1.0 1.0], pi/2)
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

julia> hflow(0.0, MMatrix{2,2}(1.0, 3.0, 2.0, 4.0), MMatrix{2,2}(0.0, 1.0, 0.0, 1.0), pi/2)   # mutable
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

julia> hflow(0.0, SMatrix{2,2}(1.0, 3.0, 2.0, 4.0), SMatrix{2,2}(0.0, 1.0, 0.0, 1.0), pi/2)   # immutable
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

Complex states ​

Complex initial conditions work with the same real Hamiltonian vector field — the scalar case:

julia
julia> hflow(0.0, 1.0 + 2.0im, 0.0 + 0.0im, pi/2)   # ≈ (0+0im, -1-2im)
(-9.676404246107521e-11 - 1.9352808492215042e-10im, -0.9999999997977443 - 1.9999999995954887im)

vectors:

julia
julia> hflow(0.0, [1.0 + 2.0im, 0.0 + 0.0im], [0.0 + 0.0im, 1.0 + 1.0im], pi/2)
(ComplexF64[-1.1098569197652709e-10 - 2.2197138395305417e-10im, 0.9999999997731949 + 0.9999999997731949im], ComplexF64[-0.9999999997731949 - 1.9999999995463897im, -1.1098569197652709e-10 - 1.1098569197652709e-10im])

and matrices:

julia
julia> hflow(0.0, [1.0+2.0im 5.0+6.0im; 3.0+4.0im 7.0+8.0im], [0.0+0.0im 1.0+1.0im; 2.0+2.0im 3.0+3.0im], pi/2)
(ComplexF64[-5.782075318723213e-11 - 1.1564150637446425e-10im 0.9999999995798117 + 0.999999999521992im; 1.999999999564368 + 1.9999999995065474im 2.9999999992019992 + 2.9999999991441815im], ComplexF64[-0.9999999998689151 - 1.9999999997378302im -4.9999999994023945 - 5.999999999271313im; -2.9999999997223874 - 3.9999999995913007im -6.999999999255864 - 7.9999999991247845im])

Automatic differentiation (ForwardDiff.Dual) ​

A ForwardDiff.Dual initial condition propagates a sensitivity through the integration — the basis for differentiating a flow with respect to its initial state/costate:

julia
julia> x0 = ForwardDiff.Dual(1.0, 1.0)    # value 1.0, seed 1.0
Dual{Nothing}(1.0,1.0)

julia> p0 = ForwardDiff.Dual(0.0, 0.0)
Dual{Nothing}(0.0,0.0)

julia> xf, pf = hflow(0.0, x0, p0, pi/2)
(Dual{Nothing}(-1.5541996200011742e-10,-1.5541996200011742e-10), Dual{Nothing}(-0.9999999997009266,-0.9999999997009266))

julia> (ForwardDiff.value(xf), ForwardDiff.value(pf))   # ≈ (0.0, -1.0)
(-1.5541996200011742e-10, -0.9999999997009266)

Vector duals work the same way:

julia
julia> x0 = [ForwardDiff.Dual(1.0, 1.0), ForwardDiff.Dual(0.0, 0.0)]
2-element Vector{ForwardDiff.Dual{Nothing, Float64, 1}}:
 Dual{Nothing}(1.0,1.0)
 Dual{Nothing}(0.0,0.0)

julia> p0 = [ForwardDiff.Dual(0.0, 0.0), ForwardDiff.Dual(1.0, 0.0)]
2-element Vector{ForwardDiff.Dual{Nothing, Float64, 1}}:
 Dual{Nothing}(0.0,0.0)
 Dual{Nothing}(1.0,0.0)

julia> hflow(0.0, x0, p0, pi/2)
(ForwardDiff.Dual{Nothing, Float64, 1}[Dual{Nothing}(-1.5541996200011742e-10,-1.5541996200011742e-10), Dual{Nothing}(0.9999999997009266,0.0)], ForwardDiff.Dual{Nothing, Float64, 1}[Dual{Nothing}(-0.9999999997009266,-0.9999999997009266), Dual{Nothing}(-1.5541996200011742e-10,0.0)])

Grid invariance

As for Flow(VectorField), the SciML integrator's internalnorm ignores dual parts, so the integration grid is identical whether the initial condition is a Real or a Dual.

Recommended: differentiate the flow, not one hand-seeded Dual ​

As for Flow(VectorField), wrap the flow call in an outer ForwardDiff.jacobian / ForwardDiff.gradient to get every partial in one call, rather than hand-seeding a Dual per component. Flow(HamiltonianVectorField) has no internal AD backend, so there is no nesting concern (unlike Flow(Hamiltonian)):

julia
julia> shoot(z) = collect(hflow(0.0, z[1], z[2], pi/2))
shoot (generic function with 1 method)

julia> ForwardDiff.jacobian(shoot, [1.0, 0.0])   # ≈ [0 1; -1 0]
2×2 Matrix{Float64}:
 -1.5542e-10   1.0
 -1.0         -1.5542e-10

This is the same Jacobian of the same dynamics computed by Flow(Hamiltonian) — a good cross-check that both constructors agree.

In-place Hamiltonian vector fields ​

An in-place Hamiltonian vector field (dx, dp, x, p) -> (dx .= p; dp .= -x) produces the same results. The state and costate are concatenated via vcat(x0, p0) into a single ODE state; the flow picks the integration path from the mutability of that concatenation. For a mutable container (Vector, MVector, Matrix, MMatrix) — and, notably, for a scalar (x0, p0) too, since vcat of two numbers always produces a mutable Vector — it is integrated truly in place:

julia
julia> hflow_ip(0.0, 1.0, 0.0, pi/2)                            # scalar: vcat(x0,p0) is mutable
(-1.5541996200011742e-10, -0.9999999997009266)

julia> hflow_ip(0.0, [1.0, 0.0], [0.0, 1.0], pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

julia> hflow_ip(0.0, MVector{2}(1.0, 0.0), MVector{2}(0.0, 1.0), pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

Immutable initial conditions (SVector / SMatrix) ​

An immutable static array cannot be written in place, so an in-place Hamiltonian vector field falls back to an out-of-place finalize RHS and emits a performance warning (note (a)). The result is still correct:

julia
julia> hflow_ip(0.0, SA[1.0, 0.0], SA[0.0, 1.0], pi/2)   # warns, then returns the correct value
┌ Warning: InPlace HamiltonianVectorField with immutable u0 (e.g. SVector): consider using an out-of-place function for better performance.
└ @ CTFlows.Systems ~/work/CTFlows.jl/CTFlows.jl/src/Systems/hamiltonian_vector_field_system.jl:341
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

For static states, prefer the out-of-place Hamiltonian vector field (hflow) shown above, or use mutable MVector / MMatrix.

See also ​