Skip to content

Commit

Permalink
Merge branch 'main' into plot
Browse files Browse the repository at this point in the history
  • Loading branch information
@ocots
ocots committed Jun 20, 2024
2 parents 424e131 + dc089a7 commit e1f9969
Show file tree
Hide file tree
Showing 5 changed files with 320 additions and 12 deletions.
142 changes: 140 additions & 2 deletions src/differential_geometry.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -515,18 +515,156 @@ julia> F0 = VectorField(x -> [x[1], x[2], (1-x[3])])
julia> F1 = VectorField(x -> [0, -x[3], x[2]])
julia> @Lie [F0, F1]([1, 2, 3])
[0, 5, 4]
#
julia> F0 = VectorField((t, x) -> [t+x[1], x[2], (1-x[3])], autonomous=false)
julia> F1 = VectorField((t, x) -> [t, -x[3], x[2]], autonomous=false)
julia> @Lie [F0, F1](1, [1, 2, 3])
#
julia> F0 = VectorField((x, v) -> [x[1]+v, x[2], (1-x[3])], variable=true)
julia> F1 = VectorField((x, v) -> [0, -x[3]-v, x[2]], variable=true)
julia> @Lie [F0, F1]([1, 2, 3], 2)
#
julia> F0 = VectorField((t, x, v) -> [t+x[1]+v, x[2], (1-x[3])], autonomous=false, variable=true)
julia> F1 = VectorField((t, x, v) -> [t, -x[3]-v, x[2]], autonomous=false, variable=true)
julia> @Lie [F0, F1](1, [1, 2, 3], 2)
#
julia> H0 = Hamiltonian((x, p) -> 0.5*(2x[1]^2+x[2]^2+p[1]^2))
julia> H1 = Hamiltonian((x, p) -> 0.5*(3x[1]^2+x[2]^2+p[2]^2))
julia> @Lie {H0, H1}([1, 2, 3], [1,0,7])
julia> @Lie {H0, H1}([1, 2, 3], [1, 0, 7])
3.0
#
julia> H0 = Hamiltonian((t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2), autonomous=false)
julia> H1 = Hamiltonian((t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2), autonomous=false)
julia> @Lie {H0, H1}(1, [1, 2, 3], [1, 0, 7])
#
julia> H0 = Hamiltonian((x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v), variable=true)
julia> H1 = Hamiltonian((x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v), variable=true)
julia> @Lie {H0, H1}([1, 2, 3], [1, 0, 7], 2)
#
julia> H0 = Hamiltonian((t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v), autonomous=false, variable=true)
julia> H1 = Hamiltonian((t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v), autonomous=false, variable=true)
julia> @Lie {H0, H1}(1, [1, 2, 3], [1, 0, 7], 2)
#
```
"""
macro Lie(expr::Expr)

fun(x) = @match (@capture(x, [a_, b_]),@capture(x, {c_, d_})) begin
(true, false) => :(Lie($a, $b))
(true, false) => :( Lie($a, $b))
(false, true) => :(Poisson($c, $d))
(false, false) => x
_ => error("internal error")
end

return esc(postwalk(fun,expr))

end

"""
$(TYPEDSIGNATURES)
Macros for Poisson brackets
# Example
```jldoctest
julia> H0 = (x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2)
julia> H1 = (x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2)
julia> @Lie {H0, H1}([1, 2, 3], [1, 0, 7]) autonomous=true variable=false
#
julia> H0 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2)
julia> H1 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2)
julia> @Lie {H0, H1}(1, [1, 2, 3], [1, 0, 7]) autonomous=false variable=false
#
julia> H0 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v)
julia> H1 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v)
julia> @Lie {H0, H1}([1, 2, 3], [1, 0, 7], 2) autonomous=true variable=true
#
julia> H0 = (t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v)
julia> H1 = (t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v)
julia> @Lie {H0, H1}(1, [1, 2, 3], [1, 0, 7], 2) autonomous=false variable=true
```
"""
macro Lie(expr::Expr, arg1, arg2)

local autonomous=true
local variable=false

@match arg1 begin
:( Autonomous ) => begin autonomous=true end
:( NonAutonomous ) => begin autonomous=false end
:( autonomous = $a ) => begin autonomous=a end
:( NonFixed ) => begin variable=true end
:( Fixed ) => begin variable=false end
:( variable = $a ) => begin variable=a end
_ => throw(IncorrectArgument("Invalid argument: " * string(arg1)))
end

@match arg2 begin
:( Autonomous ) => begin autonomous=true end
:( NonAutonomous ) => begin autonomous=false end
:( autonomous = $a ) => begin autonomous=a end
:( NonFixed ) => begin variable=true end
:( Fixed ) => begin variable=false end
:( variable = $a ) => begin variable=a end
_ => throw(IncorrectArgument("Invalid argument: " * string(arg2)))
end

fun(x) = @match (@capture(x, [a_, b_]),@capture(x, {c_, d_})) begin
#(true, false) => :( Lie($a, $b; autonomous=$autonomous, variable=$variable))
(false, true) => quote
($c isa Function && $d isa Function) ?
Poisson($c, $d; autonomous=$autonomous, variable=$variable) : Poisson($c, $d)
end
(false, false) => x
_ => error("internal error")
end

return esc(postwalk(fun,expr))

end

"""
$(TYPEDSIGNATURES)
Macros for Lie and Poisson brackets
# Example
```jldoctest
julia> H0 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2)
julia> H1 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2)
julia> @Lie {H0, H1}(1, [1, 2, 3], [1, 0, 7]) autonomous=false
#
julia> H0 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v)
julia> H1 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v)
julia> @Lie {H0, H1}([1, 2, 3], [1, 0, 7], 2) variable=true
#
```
"""
macro Lie(expr::Expr, arg)

local autonomous=true
local variable=false

@match arg begin
:( Autonomous ) => begin autonomous=true end
:( NonAutonomous ) => begin autonomous=false end
:( autonomous = $a ) => begin autonomous=a end
:( NonFixed ) => begin variable=true end
:( Fixed ) => begin variable=false end
:( variable = $a ) => begin variable=a end
_ => throw(IncorrectArgument("Invalid argument: " * string(arg)))
end

fun(x) = @match (@capture(x, [a_, b_]),@capture(x, {c_, d_})) begin
#(true, false) => :( Lie($a, $b; autonomous=$autonomous, variable=$variable))
(false, true) => quote
($c isa Function && $d isa Function) ?
Poisson($c, $d; autonomous=$autonomous, variable=$variable) : Poisson($c, $d)
end
(false, false) => x
_ => error("internal error")
end

return esc(postwalk(fun,expr))

end
2 changes: 1 addition & 1 deletion src/model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -536,7 +536,7 @@ function constraint!(
end

# bounds
(!isnothing(lb) && isnothing(ub)) && (ub = Inf*(size(lb,1) == 1 ? 1 : ones(eltype(ub), size(ub,1))))
(!isnothing(lb) && isnothing(ub)) && (ub = Inf*(size(lb,1) == 1 ? 1 : ones(eltype(lb), size(lb,1))))
( isnothing(lb) && !isnothing(ub)) && (lb = -Inf*(size(ub,1) == 1 ? 1 : ones(eltype(ub), size(ub,1))))

# range
Expand Down
14 changes: 7 additions & 7 deletions src/onepass.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ end

__init_aliases() = begin
al = OrderedDict{Symbol, Union{Real, Symbol, Expr}}()
for i ∈ 1:9 al[Symbol(:R, ctupperscripts(i))] = :( R^$i ) end
for i ∈ 1:9 al[Symbol(:R, ctupperscripts(i))] = :( R^$i ) end
al
end

Expand Down Expand Up @@ -88,15 +88,15 @@ parse!(p, ocp, e; log=false) = begin
end
# variable
:( $v ∈ R^$q, variable ) => p_variable!(p, ocp, v, q; log)
:( $v ∈ R , variable ) => p_variable!(p, ocp, v ; log)
:( $v ∈ R , variable ) => p_variable!(p, ocp, v, 1; log)
# time
:( $t ∈ [ $t0, $tf ], time ) => p_time!(p, ocp, t, t0, tf; log)
# state
:( $x ∈ R^$n, state ) => p_state!(p, ocp, x, n; log)
:( $x ∈ R , state ) => p_state!(p, ocp, x ; log)
:( $x ∈ R , state ) => p_state!(p, ocp, x, 1; log)
# control
:( $u ∈ R^$m, control ) => p_control!(p, ocp, u, m; log)
:( $u ∈ R , control ) => p_control!(p, ocp, u ; log)
:( $u ∈ R , control ) => p_control!(p, ocp, u, 1; log)
# dynamics
:( ∂($x)($t) == $e1 ) => p_dynamics!(p, ocp, x, t, e1 ; log)
:( ∂($x)($t) == $e1, $label ) => p_dynamics!(p, ocp, x, t, e1, label; log)
Expand Down Expand Up @@ -163,7 +163,7 @@ parse!(p, ocp, e; log=false) = begin
end
end

p_variable!(p, ocp, v, q=1; components_names=nothing, log=false) = begin
p_variable!(p, ocp, v, q; components_names=nothing, log=false) = begin
log && println("variable: $v, dim: $q")
v isa Symbol || return __throw("forbidden variable name: $v", p.lnum, p.line)
p.v = v
Expand Down Expand Up @@ -224,7 +224,7 @@ p_time!(p, ocp, t, t0, tf; log=false) = begin
__wrap(code, p.lnum, p.line)
end

p_state!(p, ocp, x, n=1; components_names=nothing, log=false) = begin
p_state!(p, ocp, x, n; components_names=nothing, log=false) = begin
log && println("state: $x, dim: $n")
x isa Symbol || return __throw("forbidden state name: $x", p.lnum, p.line)
p.x = x
Expand All @@ -242,7 +242,7 @@ p_state!(p, ocp, x, n=1; components_names=nothing, log=false) = begin
end
end

p_control!(p, ocp, u, m=1; components_names=nothing, log=false) = begin
p_control!(p, ocp, u, m; components_names=nothing, log=false) = begin
log && println("control: $u, dim: $m")
u isa Symbol || return __throw("forbidden control name: $u", p.lnum, p.line)
p.u = u
Expand Down
65 changes: 65 additions & 0 deletions test/test_differential_geometry.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -559,6 +559,71 @@ function test_differential_geometry()
Test.@test P011_(t, x, p, v) ≈ P011__(t, x, p, v) atol=1e-6
end

@testset "poisson macro with functions" begin
# parameters
t = 1
x = [1, 2, 3]
p = [1,0,7]
Γ = 2
γ = 1
δ = γ-Γ
v = 2

# autonomous
H0 = (x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2)
H1 = (x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2)
P01 = Poisson(H0, H1)
P011 = Poisson(P01, H1)
P01_= @Lie {H0, H1}
P011_= @Lie {{H0, H1}, H1}
Test.@test P01(x, p) ≈ P01_(x, p) atol=1e-6
Test.@test P011(x, p) ≈ P011_(x, p) atol=1e-6
get_H0 = () -> H0
P011__ = @Lie {{get_H0(), H1}, H1}
Test.@test P011_(x, p) ≈ P011__(x, p) atol=1e-6

# nonautonomous
H0 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[1]^2)
H1 = (t, x, p) -> 0.5*(x[1]^2+x[2]^2+p[2]^2)
P01 = Poisson(H0, H1; autonomous=false)
P011 = Poisson(P01, H1)
P01_val = @Lie {H0, H1}(t, x, p) autonomous=false
P01_= @Lie {H0, H1} autonomous=false
P011_= @Lie {{H0, H1}, H1} autonomous=false
Test.@test P01(t, x, p) ≈ P01_(t, x, p) atol=1e-6
Test.@test P01_val ≈ P01_(t, x, p) atol=1e-6
Test.@test P011(t, x, p) ≈ P011_(t, x, p) atol=1e-6
get_H0 = () -> H0
P011__ = @Lie {{get_H0(), H1}, H1} NonAutonomous
Test.@test P011_(t, x, p) ≈ P011__(t, x, p) atol=1e-6

# autonomous nonfixed
H0 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v)
H1 = (x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v)
P01 = Poisson(H0, H1; variable=true)
P011 = Poisson(P01, H1)
P01_= @Lie {H0, H1} variable=true
P011_= @Lie {{H0, H1}, H1} variable=true
Test.@test P01(x, p, v) ≈ P01_(x, p, v) atol=1e-6
Test.@test P011(x, p, v) ≈ P011_(x, p, v) atol=1e-6
get_H0 = () -> H0
P011__ = @Lie {{get_H0(), H1}, H1} NonFixed
Test.@test P011_(x, p, v) ≈ P011__(x, p, v) atol=1e-6

# nonautonomous nonfixed
H0 = (t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[1]^2+v)
H1 = (t, x, p, v) -> 0.5*(x[1]^2+x[2]^2+p[2]^2+v)
P01 = Poisson(H0, H1; autonomous=false, variable=true)
P011 = Poisson(P01, H1)
P01_= @Lie {H0, H1} autonomous=false variable=true
P011_= @Lie {{H0, H1}, H1} NonAutonomous NonFixed
Test.@test P01(t, x, p, v) ≈ P01_(t, x, p,v ) atol=1e-6
Test.@test P011(t, x, p, v) ≈ P011_(t, x, p, v) atol=1e-6
get_H0 = () -> H0
P011__ = @Lie {{get_H0(), H1}, H1} NonAutonomous NonFixed
Test.@test P011_(t, x, p, v) ≈ P011__(t, x, p, v) atol=1e-6
end

@testset "lie and poisson macros operation (sum,diff,...)" begin
# parameters
t = 1
Expand Down
Loading

0 comments on commit e1f9969

Please sign in to comment.