Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

foo #20

Merged
merged 1 commit into from
Jun 22, 2024
Merged

foo #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
6 changes: 3 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,16 +9,16 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"

[weakdeps]
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[extensions]
CTFlowsODE = "DifferentialEquations"
CTFlowsODE = "OrdinaryDiffEq"
CTFlowsPlots = "Plots"

[compat]
CTBase = "0.10"
DifferentialEquations = "7.13"
OrdinaryDiffEq = "6.84"
DocStringExtensions = "0.9"
Plots = "1.38"
julia = "1.8"
2 changes: 1 addition & 1 deletion docs/src/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,4 @@ The `CTFlows.jl` package is part of the [control-toolbox ecosystem](https://gith

It aims to provide tools to solve [mathematical flows](https://en.wikipedia.org/w/index.php?title=Flow_(mathematics)&oldid=1147546136#Flows_of_vector_fields_on_manifolds) of vector fields, and in particular [Hamiltonian vector fields](https://en.wikipedia.org/w/index.php?title=Hamiltonian_vector_field&oldid=1065470192) directly from the definition of the Hamiltonian, using automatic differentiation to construct the assiocated Hamiltonian vector field.

The flow is then computed thanks to [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package.
The flow is then computed thanks to [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) package.
2 changes: 1 addition & 1 deletion ext/CTFlowsODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ module CTFlowsODE

using CTBase
using CTFlows
using DifferentialEquations
using OrdinaryDiffEq
using DocStringExtensions
using MLStyle
#
Expand Down
6 changes: 3 additions & 3 deletions ext/function.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""
$(TYPEDSIGNATURES)

Returns a function that solves any ODE problem with DifferentialEquations.
Returns a function that solves any ODE problem with OrdinaryDiffEq.
"""
function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)

Expand All @@ -10,13 +10,13 @@ function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)

# ode
ode = isnothing(v) ? DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan) : DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan, v)
ode = isnothing(v) ? OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan) : OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan, v)

# jumps and callbacks
cb, t_stops_all = __callbacks(callback, jumps, nothing, _t_stops_interne, tstops)

# solve
sol = DifferentialEquations.solve(ode,
sol = OrdinaryDiffEq.solve(ode,
alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb;
kwargs_Flow..., kwargs...)

Expand Down
4 changes: 2 additions & 2 deletions ext/hamiltonian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,14 +10,14 @@ function hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)

# ode
ode = DifferentialEquations.ODEProblem(DiffEqRHS, [x0; p0], tspan, v)
ode = OrdinaryDiffEq.ODEProblem(DiffEqRHS, [x0; p0], tspan, v)

# jumps and callbacks
n = size(x0, 1)
cb, t_stops_all = __callbacks(callback, jumps, rg(n+1, 2n), _t_stops_interne, tstops)

# solve
sol = DifferentialEquations.solve(ode,
sol = OrdinaryDiffEq.solve(ode,
alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb;
kwargs_Flow..., kwargs...)

Expand Down
4 changes: 2 additions & 2 deletions ext/vector_field.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@ function vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
jumps, _t_stops_interne, DiffEqRHS, tstops=__tstops(), callback=__callback(), kwargs...)

# ode
ode = DifferentialEquations.ODEProblem(DiffEqRHS, x0, tspan, v)
ode = OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan, v)

# jumps and callbacks
cb, t_stops_all = __callbacks(callback, jumps, nothing, _t_stops_interne, tstops)

# solve
sol = DifferentialEquations.solve(ode,
sol = OrdinaryDiffEq.solve(ode,
alg=alg, abstol=abstol, reltol=reltol, saveat=saveat, tstops=t_stops_all, callback=cb;
kwargs_Flow..., kwargs...)

Expand Down
96 changes: 48 additions & 48 deletions test/test_concatenation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,38 +45,38 @@ function test_concatenation()
# one flow is used because t1 > tf
f = f1 * (2tf, f2)
xf, pf = f(t0, x0, p0, tf)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# two flows: going back
f = f1 * ((t0+tf)/2, f2)
xf, pf = f(t0, x0, p0, tf)
@test xf ≈ x0 atol = 1e-5
@test pf ≈ p0 atol = 1e-5
Test.@test xf ≈ x0 atol = 1e-5
Test.@test pf ≈ p0 atol = 1e-5

# three flows: go forward
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# autonomous and nonautonomous
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# on a grid
f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
N = 100; saveat = range(t0, tf, N)
sol = f((t0, tf), x0, p0, saveat=saveat)
xf = sol.u[end][1:n]
pf = sol.u[end][n+1:2n]
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
zspan = sol.u
zspan_sol = z_sol.(sol.t)
@test zspan ≈ zspan_sol atol = 1e-5
Test.@test zspan ≈ zspan_sol atol = 1e-5
end

@testset "Hamiltonian vector field" begin
Expand All @@ -89,38 +89,38 @@ function test_concatenation()
# one flow is used because t1 > tf
f = f1 * (2tf, f2)
xf, pf = f(t0, x0, p0, tf)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# two flows: going back
f = f1 * ((t0+tf)/2, f2)
xf, pf = f(t0, x0, p0, tf)
@test xf ≈ x0 atol = 1e-5
@test pf ≈ p0 atol = 1e-5
Test.@test xf ≈ x0 atol = 1e-5
Test.@test pf ≈ p0 atol = 1e-5

# three flows: go forward
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# autonomous and nonautonomous
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# on a grid
f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
N = 100; saveat = range(t0, tf, N)
sol = f((t0, tf), x0, p0, saveat=saveat)
xf = sol.u[end][1:n]
pf = sol.u[end][n+1:2n]
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
zspan = sol.u
zspan_sol = z_sol.(sol.t)
@test zspan ≈ zspan_sol atol = 1e-5
Test.@test zspan ≈ zspan_sol atol = 1e-5
end

@testset "Vector field" begin
Expand All @@ -133,34 +133,34 @@ function test_concatenation()
# one flow is used because t1 > tf
f = f1 * (2tf, f2)
zf = f(t0, [x0; p0], tf)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# two flows: going back
f = f1 * ((t0+tf)/2, f2)
zf = f(t0, [x0; p0], tf)
@test zf ≈ [x0; p0] atol = 1e-5
Test.@test zf ≈ [x0; p0] atol = 1e-5

# three flows: go forward
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
zf = f(t0, [x0; p0], tf+(t0+tf)/2)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# autonomous and nonautonomous
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
zf = f(t0, [x0; p0], tf+(t0+tf)/2)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# on a grid
f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
N = 100; saveat = range(t0, tf, N)
sol = f((t0, tf), [x0; p0], saveat=saveat)
xf = sol.u[end][1:n]
pf = sol.u[end][n+1:2n]
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
zspan = sol.u
zspan_sol = z_sol.(sol.t)
@test zspan ≈ zspan_sol atol = 1e-5
Test.@test zspan ≈ zspan_sol atol = 1e-5

end

Expand All @@ -174,34 +174,34 @@ function test_concatenation()
# one flow is used because t1 > tf
f = f1 * (2tf, f2)
zf = f(t0, [x0; p0], tf)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# two flows: going back
f = f1 * ((t0+tf)/2, f2)
zf = f(t0, [x0; p0], tf)
@test zf ≈ [x0; p0] atol = 1e-5
Test.@test zf ≈ [x0; p0] atol = 1e-5

# three flows: go forward
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f1)
zf = f(t0, [x0; p0], tf+(t0+tf)/2)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# autonomous and nonautonomous
f = f1 * ((t0+tf)/4, f2) * ((t0+tf)/2, f3)
zf = f(t0, [x0; p0], tf+(t0+tf)/2)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

# on a grid
f = f1 * ((t0+tf)/4, f1) * ((t0+tf)/2, f1)
N = 100; saveat = range(t0, tf, N)
sol = f((t0, tf), [x0; p0], saveat=saveat)
xf = sol.u[end][1:n]
pf = sol.u[end][n+1:2n]
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
zspan = sol.u
zspan_sol = z_sol.(sol.t)
@test zspan ≈ zspan_sol atol = 1e-5
Test.@test zspan ≈ zspan_sol atol = 1e-5

end

Expand All @@ -213,16 +213,16 @@ function test_concatenation()
f3 = Flow(Hamiltonian(H3, autonomous=false))
f = f1 * ((t0+tf)/4, [0, 0], f2) * ((t0+tf)/2, f3)
xf, pf = f(t0, x0, p0, tf+(t0+tf)/2)
@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
Test.@test pf ≈ [p1_sol(tf), p2_sol(tf)] atol = 1e-5

# vector field
f1 = Flow(VectorField(V1))
f2 = Flow(VectorField(V2))
f3 = Flow(VectorField(V3, autonomous=false))
f = f1 * ((t0+tf)/4, [0, 0, 0, 0], f2) * ((t0+tf)/2, f3)
zf = f(t0, [x0; p0], tf+(t0+tf)/2)
@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5

end

Expand Down Expand Up @@ -257,22 +257,22 @@ function test_concatenation()
# test
N = 100
tspan = range(0, 10, N)
@test norm([sol(t)[1]-sol2(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
@test norm([sol(t)[1]-sol3(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
Test.@test norm([sol(t)[1]-sol2(t) for t ∈ tspan])/N ≈ 0 atol=1e-3
Test.@test norm([sol(t)[1]-sol3(t) for t ∈ tspan])/N ≈ 0 atol=1e-3

# -------
f = Flow(Hamiltonian((x, p) -> 0.5p^2))
fc = f * (1, 1, f) * (1.5, f) * (2, 1, f) * (2.5, f) * (3, 1, f) * (3.5, f) * (4, 1, f)
xf, pf = fc(0, 0, 0, 5)
@test xf ≈ 10 atol=1e-6
@test pf ≈ 4 atol=1e-6
Test.@test xf ≈ 10 atol=1e-6
Test.@test pf ≈ 4 atol=1e-6

# -------
f = Flow(HamiltonianVectorField((x, p) -> ([p[1], 0], [0, 0])))
fc = f * (1, [1, 0], f) * (1.5, f) * (2, [1, 0], f) * (2.5, f) * (3, [1, 0], f) * (3.5, f) * (4, [1, 0], f)
xf, pf = fc(0, [0, 0], [0, 0], 5)
@test xf[1] ≈ 10 atol=1e-6
@test pf[1] ≈ 4 atol=1e-6
Test.@test xf[1] ≈ 10 atol=1e-6
Test.@test pf[1] ≈ 4 atol=1e-6

end

Expand All @@ -287,8 +287,8 @@ function test_concatenation()
f = Flow(ocp, (x, p) -> [p[1]/2, 0])
fc = f * (1, [1, 0], f) * (1.5, f) * (2, [1, 0], f) * (2.5, f) * (3, [1, 0], f) * (3.5, f) * (4, [1, 0], f)
xf, pf = fc(0, [0, 0], [0, 0], 5)
@test xf[1] ≈ 10 atol=1e-6
@test pf[1] ≈ 4 atol=1e-6
Test.@test xf[1] ≈ 10 atol=1e-6
Test.@test pf[1] ≈ 4 atol=1e-6
end

@testset "Concat OCP" begin
Expand All @@ -307,7 +307,7 @@ function test_concatenation()
t1 = -log(p0)
f = f0 * (t1, f1)
xf_, pf = f(0, -1, p0, 1)
@test xf_ ≈ 0 atol=1e-6
Test.@test xf_ ≈ 0 atol=1e-6
end

end
18 changes: 9 additions & 9 deletions test/test_default.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
function test_default()

@test CTFlows.__abstol() isa Real
@test CTFlows.__abstol() > 0
@test CTFlows.__abstol() < 1
@test CTFlows.__reltol() isa Real
@test CTFlows.__reltol() > 0
@test CTFlows.__reltol() < 1
@test CTFlows.__saveat() isa Vector
@test isnothing(CTFlows.__alg())
@test CTFlows.__tstops() isa Vector{<:Time}
Test.@test CTFlows.__abstol() isa Real
Test.@test CTFlows.__abstol() > 0
Test.@test CTFlows.__abstol() < 1
Test.@test CTFlows.__reltol() isa Real
Test.@test CTFlows.__reltol() > 0
Test.@test CTFlows.__reltol() < 1
Test.@test CTFlows.__saveat() isa Vector
Test.@test isnothing(CTFlows.__alg())
Test.@test CTFlows.__tstops() isa Vector{<:Time}

end
Loading
Loading