diff --git a/Project.toml b/Project.toml index c509ca8..7592eb0 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,10 @@ authors = ["Olivier Cots "] version = "0.8.9" [deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd" CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d" +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -19,12 +21,14 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" CTFlowsODE = "OrdinaryDiffEq" [compat] +ADTypes = "1" CTBase = "0.16" CTModels = "0.6" +DifferentiationInterface = "0.7" DocStringExtensions = "0.9" ForwardDiff = "0.10, 1.0" LinearAlgebra = "1" MLStyle = "0.4" MacroTools = "0.5" -OrdinaryDiffEq = "6" +OrdinaryDiffEq = "6, 7" julia = "1.10" diff --git a/src/CTFlows.jl b/src/CTFlows.jl index b316c76..2c88e20 100644 --- a/src/CTFlows.jl +++ b/src/CTFlows.jl @@ -7,6 +7,9 @@ using DocStringExtensions using MLStyle using MacroTools: @capture, postwalk, striplines using ForwardDiff: ForwardDiff +using ADTypes: AbstractADType +using DifferentiationInterface +using DifferentiationInterface: AutoForwardDiff, derivative, gradient # to be extended Flow(args...; kwargs...) = throw(CTBase.ExtensionError(:OrdinaryDiffEq)) @@ -18,11 +21,4 @@ include("utils.jl") include("differential_geometry.jl") include("optimal_control_problem_utils.jl") -# -# export VectorField -# export Hamiltonian -# export HamiltonianLift -# export HamiltonianVectorField -# export Flow - end diff --git a/src/default.jl b/src/default.jl index b21d2b0..2d1ba3d 100644 --- a/src/default.jl +++ b/src/default.jl @@ -44,3 +44,22 @@ Used to check whether the problem is parameterized by an external vector `v`. function __variable(ocp::CTModels.Model)::Bool return CTModels.variable_dimension(ocp) > 0 end + +""" +$(TYPEDSIGNATURES) + +Return the default automatic differentiation backend. + +Used for computing derivatives in differential geometry operations (Lie derivatives, Lie brackets, etc.). + +# Returns +- An instance of `AutoForwardDiff()` from DifferentiationInterface.jl. + +# Example +```julia-repl +julia> backend = __backend() +julia> typeof(backend) +AutoForwardDiff +``` +""" +__backend() = AutoForwardDiff() diff --git a/src/differential_geometry.jl b/src/differential_geometry.jl index 8b53564..06e1689 100644 --- a/src/differential_geometry.jl +++ b/src/differential_geometry.jl @@ -1,671 +1,912 @@ +# ============================================================================== +# CTFlows.jl - Differential Geometry Module (V3) +# ============================================================================== +# This file contains differential geometry tools for control theory. +# V3 refactoring: wrapper-free API using pure Functions + +# ============================================================================== +# Prefix System for Differential Geometry Functions +# ============================================================================== + +""" +Prefix reference for differential geometry functions in the `@Lie` macro. + +This constant stores a Symbol indicating which module should be used when the `@Lie` macro +expands Lie bracket `[X, Y]` and Poisson bracket `{H, G}` expressions. By default, it points +to `:CTFlows`, but can be changed via [`diffgeo_prefix!`](@ref) to support custom modules. + +# Example + +```julia-repl +julia> using CTFlows + +julia> CTFlows.diffgeo_prefix() +:CTFlows + +julia> CTFlows.diffgeo_prefix!(:MyModule) + +julia> CTFlows.diffgeo_prefix() +:MyModule +``` + +See also: [`diffgeo_prefix`](@ref), [`diffgeo_prefix!`](@ref), [`@Lie`](@ref) +""" +const DIFFGEO_PREFIX = Ref(:CTFlows) + +""" +$(TYPEDSIGNATURES) + +Get the current differential geometry module prefix used by the `@Lie` macro. + +# Returns + +- `Symbol`: The current prefix (default: `:CTFlows`) + +# Example + +```julia-repl +julia> using CTFlows + +julia> CTFlows.diffgeo_prefix() +:CTFlows +``` + +See also: [`diffgeo_prefix!`](@ref), [`@Lie`](@ref) +""" +diffgeo_prefix() = DIFFGEO_PREFIX[] + """ $(TYPEDSIGNATURES) -Construct the Hamiltonian lift of a `VectorField`. +Set the differential geometry module prefix for the `@Lie` macro. + +This changes which module the `@Lie` macro will reference when expanding bracket expressions. +Useful when integrating with custom differential geometry implementations. # Arguments -- `X::VectorField`: The vector field to lift. Its signature determines if it is autonomous and/or variable. -# Returns -- A `HamiltonianLift` callable object representing the Hamiltonian lift of `X`. +- `p::Symbol`: The new module prefix (e.g., `:MyCustomModule`) + +# Example -# Examples ```julia-repl -julia> HL = Lift(VectorField(x -> [x[1]^2, x[2]^2], autonomous=true, variable=false)) -julia> HL([1, 0], [0, 1]) # returns 0 +julia> using CTFlows -julia> HL2 = Lift(VectorField((t, x, v) -> [t + x[1]^2, x[2]^2 + v], autonomous=false, variable=true)) -julia> HL2(1, [1, 0], [0, 1], 1) # returns 1 +julia> CTFlows.diffgeo_prefix!(:MyModule) -julia> H = Lift(x -> 2x) -julia> H(1, 1) # returns 2 +julia> CTFlows.diffgeo_prefix() +:MyModule -julia> H2 = Lift((t, x, v) -> 2x + t - v, autonomous=false, variable=true) -julia> H2(1, 1, 1, 1) # returns 2 +julia> # Now @Lie macros will expand to MyModule.ad(...) instead of CTFlows.ad(...) -# Alternative syntax using symbols for autonomy and variability -julia> H3 = Lift((t, x, v) -> 2x + t - v, NonAutonomous, NonFixed) -julia> H3(1, 1, 1, 1) # returns 2 +julia> CTFlows.diffgeo_prefix!(:CTFlows) # Reset to default ``` + +See also: [`diffgeo_prefix`](@ref), [`@Lie`](@ref) """ -function Lift(X::VectorField)::HamiltonianLift - return HamiltonianLift(X) +function diffgeo_prefix!(p::Symbol) + DIFFGEO_PREFIX[] = p + return nothing end +# ============================================================================== +# Unified ad() Function - Lie Derivative and Lie Bracket +# ============================================================================== + """ $(TYPEDSIGNATURES) -Construct the Hamiltonian lift of a function. +Unified function for Lie derivative and Lie bracket using directional derivatives. # Arguments -- `X::Function`: The function representing the vector field. -- `autonomous::Bool=true`: Whether the function is autonomous (time-independent). -- `variable::Bool=false`: Whether the function depends on an additional variable argument. +- `X::Function`: Vector field +- `foo::Function`: Either a scalar function (for Lie derivative) or vector field (for Lie bracket) +- `backend`: Automatic differentiation backend (default: `__backend()` = `AutoForwardDiff()`) +- `autonomous::Bool`: Whether functions are time-independent (default: `__autonomous()` = `true`) +- `variable::Bool`: Whether functions depend on extra variable (default: `__variable()` = `false`) # Returns -- A callable function computing the Hamiltonian lift, -(and variants depending on `autonomous` and `variable`). +- A function computing either: + - Lie derivative `ad(X, f)(x) = ∇f(x)' * X(x)` if `foo` is scalar + - Lie bracket `ad(X, Y)(x) = J_Y(x)*X(x) - J_X(x)*Y(x)` if `foo` is vector -# Details -Depending on the `autonomous` and `variable` flags, the returned function has one of the following call signatures: -- `(x, p)` if `autonomous=true` and `variable=false` -- `(x, p, v)` if `autonomous=true` and `variable=true` -- `(t, x, p)` if `autonomous=false` and `variable=false` -- `(t, x, p, v)` if `autonomous=false` and `variable=true` +# Mathematical Approach +Uses directional derivatives: `D_X foo(x) = d/dt [foo(x + t*X(x))]|_{t=0}` # Examples ```julia-repl -julia> H = Lift(x -> 2x) -julia> H(1, 1) # returns 2 +julia> using CTFlows + +# Lie derivative of a scalar function +julia> X(x) = [x[2], -x[1]] + +julia> f(x) = x[1]^2 + x[2]^2 -julia> H2 = Lift((t, x, v) -> 2x + t - v, autonomous=false, variable=true) -julia> H2(1, 1, 1, 1) # returns 2 +julia> Lf = CTFlows.ad(X, f) + +julia> Lf([1.0, 2.0]) +0.0 + +# Lie bracket of two vector fields +julia> Y(x) = [x[1], x[2]] + +julia> Z = CTFlows.ad(X, Y) + +julia> Z([1.0, 2.0]) +2-element Vector{Float64}: + -1.0 + 2.0 ``` -""" -function Lift( - X::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -)::Function - return @match (autonomous, variable) begin - (true, false) => (x, p) -> p' * X(x) - (true, true) => (x, p, v) -> p' * X(x, v) - (false, false) => (t, x, p) -> p' * X(t, x) - _ => (t, x, p, v) -> p' * X(t, x, v) - end -end -# --------------------------------------------------------------------------- -# Lie derivative of a scalar function along a vector field or a function: L_X(f) = X⋅f +See also: [`Lift`](@ref), [`Poisson`](@ref), [`@Lie`](@ref) +""" +function ad( + X::Function, + foo::Function; + backend::AbstractADType=__backend(), + autonomous::Bool=__autonomous(), + variable::Bool=__variable(), +) + TD = autonomous ? Autonomous : NonAutonomous + VD = variable ? NonFixed : Fixed + return ad(X, foo, TD, VD; backend=backend) +end -# (postulate) -# (X⋅f)(x) = f'(x)⋅X(x) -# (X⋅f)(t, x) = ∂ₓf(t, x)⋅X(t, x) -# (g⋅f)(x) = f'(x)⋅G(x) with G the vector of g """ -Lie derivative of a scalar function along a vector field in the autonomous case. +$(TYPEDSIGNATURES) + +Unified function for Lie derivative and Lie bracket with explicit type parameters. + +This method accepts types directly for better performance through compile-time dispatch. +Prefer this form when the time and variable dependence are known at compile time. + +# Arguments + +- `X::Function`: Vector field +- `foo::Function`: Either a scalar function (for Lie derivative) or vector field (for Lie bracket) +- `TD::Type{<:TimeDependence}`: Time dependence type (`Autonomous` or `NonAutonomous`) +- `VD::Type{<:VariableDependence}`: Variable dependence type (`Fixed` or `NonFixed`) +- `backend::AbstractADType`: Automatic differentiation backend (default: `__backend()`) + +# Returns + +- A function with the appropriate signature based on `TD` and `VD` + +# Example -Example: ```julia-repl -julia> φ = x -> [x[2], -x[1]] -julia> X = VectorField(φ) -julia> f = x -> x[1]^2 + x[2]^2 -julia> (X⋅f)([1, 2]) -0 +julia> using CTFlows + +julia> X(x) = [x[2], -x[1]] + +julia> f(x) = x[1]^2 + x[2]^2 + +# Type-based dispatch for better performance +julia> Lf = CTFlows.ad(X, f, CTFlows.Autonomous, CTFlows.Fixed) + +julia> Lf([1.0, 2.0]) +0.0 ``` + +See also: [`ad`](@ref) """ -function ⋅( - X::VectorField{<:Function,Autonomous,<:VariableDependence}, f::Function -)::Function - return (x, args...) -> ctgradient(y -> f(y, args...), x)' * X(x, args...) +function ad( + X::Function, + foo::Function, + ::Type{TD}, + ::Type{VD}; + backend::AbstractADType=__backend(), +) where {TD<:TimeDependence,VD<:VariableDependence} + return _ad(X, foo, backend, TD, VD) end +# ============================================================================== +# Internal ad() implementations with type dispatch +# ============================================================================== + """ -Lie derivative of a scalar function along a vector field in the nonautonomous case. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> φ = (t, x, v) -> [t + x[2] + v[1], -x[1] + v[2]] -julia> X = VectorField(φ, NonAutonomous, NonFixed) -julia> f = (t, x, v) -> t + x[1]^2 + x[2]^2 -julia> (X⋅f)(1, [1, 2], [2, 1]) -10 -``` +Internal implementation of `ad` for autonomous, time-fixed functions. + +Computes the Lie derivative or Lie bracket for functions with signature `(x) -> ...`. + +This is the most common case: time-independent vector fields and functions without +extra parameters. + +# Implementation + +Uses directional derivative: `d/dt [foo(x + t*X(x))]|_{t=0}` via automatic differentiation. + +# Returns + +- A function `(x) -> result` where `result` is: + - A scalar if `foo` is a scalar function (Lie derivative) + - A vector if `foo` is a vector field (Lie bracket) + +# Note + +This is an internal function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad_result`](@ref) """ -function ⋅( - X::VectorField{<:Function,NonAutonomous,<:VariableDependence}, f::Function -)::Function - return (t, x, args...) -> ctgradient(y -> f(t, y, args...), x)' * X(t, x, args...) +function _ad( + X::Function, foo::Function, backend, ::Type{Autonomous}, ::Type{Fixed} +) + return function (x) + X_x = X(x) + g(t) = foo(x + t * X_x) + dfoo = derivative(g, backend, 0.0) + return _ad_result(X, foo, dfoo, x, X_x, backend) + end end """ -Lie derivative of a scalar function along a function (considered autonomous and non-variable). +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> φ = x -> [x[2], -x[1]] -julia> f = x -> x[1]^2 + x[2]^2 -julia> (φ⋅f)([1, 2]) -0 -julia> φ = (t, x, v) -> [t + x[2] + v[1], -x[1] + v[2]] -julia> f = (t, x, v) -> t + x[1]^2 + x[2]^2 -julia> (φ⋅f)(1, [1, 2], [2, 1]) -MethodError -``` +Internal implementation of `ad` for autonomous, variable-dependent functions. + +Computes the Lie derivative or Lie bracket for functions with signature `(x, v) -> ...`. + +This variant is used when the vector field or scalar function depends on an extra +variable `v` (e.g., control parameter, model parameter). + +# Implementation + +Uses directional derivative: `d/dt [foo(x + t*X(x,v), v)]|_{t=0}` via automatic differentiation. + +# Returns + +- A function `(x, v) -> result` where `result` is: + - A scalar if `foo` is a scalar function (Lie derivative) + - A vector if `foo` is a vector field (Lie bracket) + +# Note + +This is an internal function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad_result`](@ref) """ -function ⋅(X::Function, f::Function)::Function - return ⋅(VectorField(X, Autonomous, Fixed), f) +function _ad( + X::Function, foo::Function, backend, ::Type{Autonomous}, ::Type{NonFixed} +) + return function (x, v) + X_x = X(x, v) + g(t) = foo(x + t * X_x, v) + dfoo = derivative(g, backend, 0.0) + return _ad_result(X, foo, dfoo, x, X_x, backend, v) + end end """ -Lie derivative of a scalar function along a vector field. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> φ = x -> [x[2], -x[1]] -julia> X = VectorField(φ) -julia> f = x -> x[1]^2 + x[2]^2 -julia> Lie(X,f)([1, 2]) -0 -julia> φ = (t, x, v) -> [t + x[2] + v[1], -x[1] + v[2]] -julia> X = VectorField(φ, NonAutonomous, NonFixed) -julia> f = (t, x, v) -> t + x[1]^2 + x[2]^2 -julia> Lie(X, f)(1, [1, 2], [2, 1]) -10 -``` +Internal implementation of `ad` for non-autonomous, time-fixed functions. + +Computes the Lie derivative or Lie bracket for functions with signature `(t, x) -> ...`. + +This variant is used when the vector field or scalar function explicitly depends on time `t`. + +# Implementation + +Uses directional derivative: `d/ds [foo(t, x + s*X(t,x))]|_{s=0}` via automatic differentiation. +Note that time `t` is held constant during the directional derivative. + +# Returns + +- A function `(t, x) -> result` where `result` is: + - A scalar if `foo` is a scalar function (Lie derivative) + - A vector if `foo` is a vector field (Lie bracket) + +# Note + +This is an internal function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad_result`](@ref) """ -Lie(X::VectorField, f::Function)::Function = X ⋅ f +function _ad( + X::Function, foo::Function, backend, ::Type{NonAutonomous}, ::Type{Fixed} +) + return function (t, x) + X_x = X(t, x) + g(s) = foo(t, x + s * X_x) + dfoo = derivative(g, backend, 0.0) + return _ad_result(X, foo, dfoo, x, X_x, backend, t) + end +end """ -Lie derivative of a scalar function along a function with specified dependencies. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> φ = x -> [x[2], -x[1]] -julia> f = x -> x[1]^2 + x[2]^2 -julia> Lie(φ,f)([1, 2]) -0 -julia> φ = (t, x, v) -> [t + x[2] + v[1], -x[1] + v[2]] -julia> f = (t, x, v) -> t + x[1]^2 + x[2]^2 -julia> Lie(φ, f, autonomous=false, variable=true)(1, [1, 2], [2, 1]) -10 -``` +Internal implementation of `ad` for non-autonomous, variable-dependent functions. + +Computes the Lie derivative or Lie bracket for functions with signature `(t, x, v) -> ...`. + +This is the most general case, handling time-dependent vector fields with extra variables. + +# Implementation + +Uses directional derivative: `d/ds [foo(t, x + s*X(t,x,v), v)]|_{s=0}` via automatic differentiation. +Time `t` and variable `v` are held constant during the directional derivative. + +# Returns + +- A function `(t, x, v) -> result` where `result` is: + - A scalar if `foo` is a scalar function (Lie derivative) + - A vector if `foo` is a vector field (Lie bracket) + +# Note + +This is an internal function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad_result`](@ref) """ -function Lie( - X::Function, f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -)::Function - return Lie(VectorField(X; autonomous=autonomous, variable=variable), f) +function _ad( + X::Function, foo::Function, backend, ::Type{NonAutonomous}, ::Type{NonFixed} +) + return function (t, x, v) + X_x = X(t, x, v) + g(s) = foo(t, x + s * X_x, v) + dfoo = derivative(g, backend, 0.0) + return _ad_result(X, foo, dfoo, x, X_x, backend, t, v) + end end +# ============================================================================== +# Result computation dispatch (Lie derivative vs Lie bracket) +# ============================================================================== + """ -Partial derivative with respect to time of a function. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> ∂ₜ((t,x) -> t*x)(0,8) -8 -``` +Compute final result for Lie derivative (scalar function case). + +When `foo` is a scalar function, the directional derivative `d/dt [foo(x + t*X(x))]|_{t=0}` +directly gives the Lie derivative `∇f(x)' * X(x)`, so we just return it. + +# Arguments + +- `X::Function`: The vector field +- `foo::Function`: The scalar function +- `dfoo::Number`: The computed directional derivative (already equals Lie derivative) +- `x`: Current point +- `X_x`: Value of `X(x)` (precomputed) +- `backend`: AD backend +- `args...`: Additional arguments (time and/or variable) + +# Returns + +- `Number`: The Lie derivative value `L_X f(x) = ∇f(x)' * X(x)` + +# Note + +This is an internal dispatch function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad`](@ref) """ -∂ₜ(f) = (t, args...) -> ctgradient(y -> f(y, args...), t) +function _ad_result(X::Function, foo::Function, dfoo::Number, x, X_x, backend, args...) + return dfoo # Already ∇f(x)' * X(x) +end """ -"Directional derivative" of a vector field in the autonomous case, -used internally for computing the Lie bracket. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> X = VectorField(x -> [x[2], -x[1]]) -julia> Y = VectorField(x -> [x[1], x[2]]) -julia> (X ⅋ Y)([1, 2]) -[2, -1] -``` +Compute final result for Lie bracket (vector field case). + +When `foo` is a vector field `Y`, the directional derivative gives `J_Y(x) * X(x)`. +We must subtract `J_X(x) * Y(x)` to get the full Lie bracket `[X, Y](x)`. + +# Arguments + +- `X::Function`: The first vector field +- `foo::Function`: The second vector field (Y) +- `dfoo::AbstractVector`: The computed `J_Y(x) * X(x)` +- `x`: Current point +- `X_x`: Value of `X(x)` (precomputed) +- `backend`: AD backend +- `args...`: Additional arguments (time and/or variable) + +# Returns + +- `AbstractVector`: The Lie bracket `[X, Y](x) = J_Y(x)*X(x) - J_X(x)*Y(x)` + +# Note + +This is an internal dispatch function. Users should call [`ad`](@ref) instead. + +See also: [`ad`](@ref), [`_ad`](@ref) """ -function ⅋( - X::VectorField{<:Function,Autonomous,V}, Y::VectorField{<:Function,Autonomous,V} -)::VectorField{<:Function,Autonomous,V} where {V<:VariableDependence} - return VectorField( - (x, args...) -> if x isa ctNumber - ctgradient(y -> Y(y, args...), x) * X(x, args...) - else - ctjacobian(y -> Y(y, args...), x) * X(x, args...) - end, Autonomous, V - ) +function _ad_result( + X::Function, foo::Function, dfoo::AbstractVector, x, X_x, backend, args... +) + # dfoo = J_Y(x) * X(x) + # Compute J_X(x) * Y(x) using directional derivative + Y_x = foo(x, args...) + h(t) = X(x + t * Y_x, args...) + dX = derivative(h, backend, 0.0) + + return dfoo - dX # J_Y(x)*X(x) - J_X(x)*Y(x) end +# ============================================================================== +# Lift Function - Hamiltonian Lift (V3: returns pure Function) +# ============================================================================== + """ -"Directional derivative" of a vector field in the nonautonomous case, -used internally for computing the Lie bracket. +$(TYPEDSIGNATURES) + +Construct the Hamiltonian lift of a vector field (pure function). + +# Arguments + +- `f::Function`: Vector field function +- `autonomous::Bool`: Whether the function is autonomous (default: `__autonomous()` = `true`) +- `variable::Bool`: Whether the function depends on an additional variable (default: `__variable()` = `false`) + +# Returns + +- A callable function computing the Hamiltonian lift `H(x, p) = p' * f(x)` + +# Examples -Example: ```julia-repl -julia> X = VectorField((t, x, v) -> [t + v[1] + v[2] + x[2], -x[1]], NonFixed, NonAutonomous) -julia> Y = VectorField((t, x, v) -> [v[1] + v[2] + x[1], x[2]], NonFixed, NonAutonomous) -julia> (X ⅋ Y)(1, [1, 2], [2, 3]) -[8, -1] +julia> using CTFlows + +julia> f(x) = [x[1]^2, x[2]^2] + +julia> H = CTFlows.Lift(f) + +julia> H([1.0, 2.0], [0.5, 0.5]) +2.5 + +julia> g(t, x) = [t*x[1], x[2]] + +julia> H2 = CTFlows.Lift(g; autonomous=false) + +julia> H2(1.0, [1.0, 2.0], [0.5, 0.5]) +1.5 ``` + +See also: [`ad`](@ref), [`Poisson`](@ref) """ -function ⅋( - X::VectorField{<:Function,NonAutonomous,V}, Y::VectorField{<:Function,NonAutonomous,V} -)::VectorField{<:Function,NonAutonomous,V} where {V<:VariableDependence} - return VectorField( - (t, x, args...) -> if x isa ctNumber - ctgradient(y -> Y(t, y, args...), x) * X(t, x, args...) - else - ctjacobian(y -> Y(t, y, args...), x) * X(t, x, args...) - end, - NonAutonomous, - V, - ) +function Lift(f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable()) + TD = autonomous ? Autonomous : NonAutonomous + VD = variable ? NonFixed : Fixed + return Lift(f, TD, VD) end -# --------------------------------------------------------------------------- -# Lie bracket of two vector fields: [X, Y] = Lie(X, Y) -# [X, Y]⋅f = (X∘Y - Y∘X)⋅f = (X⅋Y - Y⅋X)⋅f """ -Lie bracket of two vector fields in the autonomous case. +$(TYPEDSIGNATURES) + +Construct the Hamiltonian lift with explicit type parameters. + +This method uses compile-time type dispatch for better performance when the +time and variable dependence are known statically. + +# Arguments + +- `f::Function`: Vector field function +- `TD::Type{<:TimeDependence}`: Time dependence type (`Autonomous` or `NonAutonomous`) +- `VD::Type{<:VariableDependence}`: Variable dependence type (`Fixed` or `NonFixed`) + +# Returns + +- A function computing `H(x, p, ...) = p' * f(x, ...)` + +# Example -Example: ```julia-repl -julia> f = x -> [x[2], 2x[1]] -julia> g = x -> [3x[2], -x[1]] -julia> X = VectorField(f) -julia> Y = VectorField(g) -julia> Lie(X, Y)([1, 2]) -[7, -14] +julia> using CTFlows + +julia> f(x) = [x[1]^2, x[2]] + +julia> H = CTFlows.Lift(f, CTFlows.Autonomous, CTFlows.Fixed) + +julia> H([1.0, 2.0], [0.5, 0.5]) +1.5 ``` + +See also: [`Lift`](@ref) """ -function Lie( - X::VectorField{<:Function,Autonomous,V}, Y::VectorField{<:Function,Autonomous,V} -)::VectorField{<:Function,Autonomous,V} where {V<:VariableDependence} - return VectorField( - (x, args...) -> (X ⅋ Y)(x, args...) - (Y ⅋ X)(x, args...), Autonomous, V - ) +function Lift( + f::Function, + ::Type{TD}, + ::Type{VD}, +) where {TD<:TimeDependence,VD<:VariableDependence} + return _Lift(f, TD, VD) end +# Internal implementations + """ -Lie bracket of two vector fields in the nonautonomous case. +$(TYPEDSIGNATURES) -Example: -```julia-repl -julia> f = (t, x, v) -> [t + x[2] + v, -2x[1] - v] -julia> g = (t, x, v) -> [t + 3x[2] + v, -x[1] - v] -julia> X = VectorField(f, NonAutonomous, NonFixed) -julia> Y = VectorField(g, NonAutonomous, NonFixed) -julia> Lie(X, Y)(1, [1, 2], 1) -[-7, 12] -``` +Internal implementation of Hamiltonian lift for autonomous, time-fixed vector fields. + +Constructs `H(x, p) = p' * f(x)` where `f` is a vector field. + +# Note + +This is an internal function. Users should call [`Lift`](@ref) instead. + +See also: [`Lift`](@ref) """ -function Lie( - X::VectorField{<:Function,NonAutonomous,V}, Y::VectorField{<:Function,NonAutonomous,V} -)::VectorField{<:Function,NonAutonomous,V} where {V<:VariableDependence} - return VectorField( - (t, x, args...) -> (X ⅋ Y)(t, x, args...) - (Y ⅋ X)(t, x, args...), NonAutonomous, V - ) -end +_Lift(f::Function, ::Type{Autonomous}, ::Type{Fixed}) = (x, p) -> p' * f(x) -# --------------------------------------------------------------------------- -# Poisson bracket of two Hamiltonian functions: {f, g} = Poisson(f, g) -# f(z) = p⋅X(x), g(z) = p⋅Y(x), z = (x, p) => {f,g}(z) = p⋅[X,Y](x) -# {f, g}(z) = g'(z)⋅F(z), where F is the Hamiltonian vector field of f -# actually, z = (x, p) """ $(TYPEDSIGNATURES) -Poisson bracket of two Hamiltonian functions (subtype of AbstractHamiltonian). Autonomous case. +Internal implementation of Hamiltonian lift for autonomous, variable-dependent vector fields. -Returns a Hamiltonian representing the Poisson bracket `{f, g}` of two autonomous Hamiltonian functions `f` and `g`. +Constructs `H(x, p, v) = p' * f(x, v)` where `f` depends on an extra variable `v`. -# Example -```julia-repl -julia> f = (x, p) -> x[2]^2 + 2x[1]^2 + p[1]^2 -julia> g = (x, p) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] -julia> F = Hamiltonian(f) -julia> G = Hamiltonian(g) -julia> Poisson(f, g)([1, 2], [2, 1]) # -20 -julia> Poisson(f, G)([1, 2], [2, 1]) # -20 -julia> Poisson(F, g)([1, 2], [2, 1]) # -20 -``` +# Note + +This is an internal function. Users should call [`Lift`](@ref) instead. + +See also: [`Lift`](@ref) """ -function Poisson( - f::AbstractHamiltonian{Autonomous,V}, g::AbstractHamiltonian{Autonomous,V} -)::Hamiltonian{<:Function,Autonomous,V} where {V<:VariableDependence} - function fg(x, p, args...) - n = size(x, 1) - ff, gg = @match n begin - 1 => (z -> f(z[1], z[2], args...), z -> g(z[1], z[2], args...)) - _ => ( - z -> f(z[1:n], z[(n + 1):2n], args...), - z -> g(z[1:n], z[(n + 1):2n], args...), - ) - end - df = ctgradient(ff, [x; p]) - dg = ctgradient(gg, [x; p]) - return df[(n + 1):2n]' * dg[1:n] - df[1:n]' * dg[(n + 1):2n] - end - return Hamiltonian(fg, Autonomous, V) -end +_Lift(f::Function, ::Type{Autonomous}, ::Type{NonFixed}) = (x, p, v) -> p' * f(x, v) """ $(TYPEDSIGNATURES) -Poisson bracket of two Hamiltonian functions. Non-autonomous case. +Internal implementation of Hamiltonian lift for non-autonomous, time-fixed vector fields. -Returns a Hamiltonian representing `{f, g}` where `f` and `g` are time-dependent. +Constructs `H(t, x, p) = p' * f(t, x)` where `f` depends on time `t`. + +# Note + +This is an internal function. Users should call [`Lift`](@ref) instead. + +See also: [`Lift`](@ref) +""" +_Lift(f::Function, ::Type{NonAutonomous}, ::Type{Fixed}) = (t, x, p) -> p' * f(t, x) + +""" +$(TYPEDSIGNATURES) + +Internal implementation of Hamiltonian lift for non-autonomous, variable-dependent vector fields. + +Constructs `H(t, x, p, v) = p' * f(t, x, v)` where `f` depends on time and an extra variable. + +# Note + +This is an internal function. Users should call [`Lift`](@ref) instead. + +See also: [`Lift`](@ref) +""" +_Lift(f::Function, ::Type{NonAutonomous}, ::Type{NonFixed}) = + (t, x, p, v) -> p' * f(t, x, v) + +# ============================================================================== +# Poisson Bracket (V3: works with pure Functions) +# ============================================================================== + +""" +$(TYPEDSIGNATURES) + +Poisson bracket of two Hamiltonian functions (pure functions). + +# Arguments + +- `H::Function`: First Hamiltonian function +- `G::Function`: Second Hamiltonian function +- `backend::AbstractADType`: Automatic differentiation backend (default: `__backend()` = `AutoForwardDiff()`) +- `autonomous::Bool`: Whether functions are time-independent (default: `__autonomous()` = `true`) +- `variable::Bool`: Whether functions depend on extra variable (default: `__variable()` = `false`) + +# Returns + +- A function computing the Poisson bracket `{H, G}(x, p) = ∇ₚH'·∇ₓG - ∇ₓH'·∇ₚG` + +# Examples -# Example ```julia-repl -julia> f = (t, x, p, v) -> t*v[1]*x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] -julia> g = (t, x, p, v) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] + t - v[2] -julia> F = Hamiltonian(f, autonomous=false, variable=true) -julia> G = Hamiltonian(g, autonomous=false, variable=true) -julia> Poisson(F, G)(2, [1, 2], [2, 1], [4, 4]) # -76 -julia> Poisson(f, g, NonAutonomous, NonFixed)(2, [1, 2], [2, 1], [4, 4]) # -76 +julia> using CTFlows + +julia> H(x, p) = 0.5 * (p[1]^2 + p[2]^2) + +julia> G(x, p) = x[1] + +julia> PB = CTFlows.Poisson(H, G) + +julia> PB([1.0, 2.0], [3.0, 4.0]) +3.0 ``` + +See also: [`Lift`](@ref), [`ad`](@ref), [`@Lie`](@ref) """ function Poisson( - f::AbstractHamiltonian{NonAutonomous,V}, g::AbstractHamiltonian{NonAutonomous,V} -)::Hamiltonian{<:Function,NonAutonomous,V} where {V<:VariableDependence} - function fg(t, x, p, args...) - n = size(x, 1) - ff, gg = @match n begin - 1 => (z -> f(t, z[1], z[2], args...), z -> g(t, z[1], z[2], args...)) - _ => ( - z -> f(t, z[1:n], z[(n + 1):2n], args...), - z -> g(t, z[1:n], z[(n + 1):2n], args...), - ) - end - df = ctgradient(ff, [x; p]) - dg = ctgradient(gg, [x; p]) - return df[(n + 1):2n]' * dg[1:n] - df[1:n]' * dg[(n + 1):2n] - end - return Hamiltonian(fg, NonAutonomous, V) + H::Function, + G::Function; + backend::AbstractADType=__backend(), + autonomous::Bool=__autonomous(), + variable::Bool=__variable(), +) + TD = autonomous ? Autonomous : NonAutonomous + VD = variable ? NonFixed : Fixed + return Poisson(H, G, TD, VD; backend=backend) end """ $(TYPEDSIGNATURES) -Poisson bracket of two HamiltonianLift vector fields. +Poisson bracket with explicit type parameters for compile-time dispatch. -Returns the HamiltonianLift corresponding to the Lie bracket of vector fields `f.X` and `g.X`. +This method uses type parameters for better performance when time and variable +dependence are known statically. + +# Arguments + +- `H::Function`: First Hamiltonian function +- `G::Function`: Second Hamiltonian function +- `TD::Type{<:TimeDependence}`: Time dependence type (`Autonomous` or `NonAutonomous`) +- `VD::Type{<:VariableDependence}`: Variable dependence type (`Fixed` or `NonFixed`) +- `backend::AbstractADType`: Automatic differentiation backend (default: `__backend()`) + +# Returns + +- A function computing `{H, G}(x, p, ...) = ∇ₚ H · ∇ₓ G - ∇ₓ H · ∇ₚ G` # Example + ```julia-repl -julia> f = x -> [x[1]^2 + x[2]^2, 2x[1]^2] -julia> g = x -> [3x[2]^2, x[2] - x[1]^2] -julia> F = Lift(f) -julia> G = Lift(g) -julia> Poisson(F, G)([1, 2], [2, 1]) # -64 - -julia> f = (t, x, v) -> [t*v[1]*x[2]^2, 2x[1]^2 + v[2]] -julia> g = (t, x, v) -> [3x[2]^2 - x[1]^2, t - v[2]] -julia> F = Lift(f, NonAutonomous, NonFixed) -julia> G = Lift(g, NonAutonomous, NonFixed) -julia> Poisson(F, G)(2, [1, 2], [2, 1], [4, 4]) # 100 +julia> using CTFlows + +julia> H(x, p) = 0.5 * (p[1]^2 + p[2]^2) + +julia> G(x, p) = x[1] + +julia> PB = CTFlows.Poisson(H, G, CTFlows.Autonomous, CTFlows.Fixed) + +julia> PB([1.0, 2.0], [3.0, 4.0]) +3.0 ``` + +See also: [`Poisson`](@ref) """ function Poisson( - f::HamiltonianLift{T,V}, g::HamiltonianLift{T,V} -)::HamiltonianLift{<:VectorField,T,V} where {T<:TimeDependence,V<:VariableDependence} - return HamiltonianLift(Lie(f.X, g.X)) + H::Function, + G::Function, + ::Type{TD}, + ::Type{VD}; + backend::AbstractADType=__backend(), +) where {TD<:TimeDependence,VD<:VariableDependence} + return _Poisson(H, G, backend, TD, VD) end +# Internal implementations + """ $(TYPEDSIGNATURES) -Poisson bracket of two functions. The time and variable dependence are specified with keyword arguments. +Internal implementation of Poisson bracket for autonomous, time-fixed Hamiltonians. -Returns a Hamiltonian computed from the functions promoted as Hamiltonians. +Computes `{H, G}(x, p) = ∇ₚ H(x,p) · ∇ₓ G(x,p) - ∇ₓ H(x,p) · ∇ₚ G(x,p)`. -# Example -```julia-repl -julia> f = (x, p) -> x[2]^2 + 2x[1]^2 + p[1]^2 -julia> g = (x, p) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] -julia> Poisson(f, g)([1, 2], [2, 1]) # -20 +# Implementation -julia> f = (t, x, p, v) -> t*v[1]*x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] -julia> g = (t, x, p, v) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] + t - v[2] -julia> Poisson(f, g, autonomous=false, variable=true)(2, [1, 2], [2, 1], [4, 4]) # -76 -``` +Uses `ctgradient` to compute partial derivatives with respect to `x` and `p` separately. + +# Note + +This is an internal function. Users should call [`Poisson`](@ref) instead. + +See also: [`Poisson`](@ref), [`ctgradient`](@ref) """ -function Poisson( - f::Function, g::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -)::Hamiltonian - return Poisson( - Hamiltonian(f; autonomous=autonomous, variable=variable), - Hamiltonian(g; autonomous=autonomous, variable=variable), - ) +function _Poisson(H::Function, G::Function, backend, ::Type{Autonomous}, ::Type{Fixed}) + return function (x, p) + grad_x_H = ctgradient(y -> H(y, p), backend, x) + grad_p_H = ctgradient(q -> H(x, q), backend, p) + grad_x_G = ctgradient(y -> G(y, p), backend, x) + grad_p_G = ctgradient(q -> G(x, q), backend, p) + return grad_p_H' * grad_x_G - grad_x_H' * grad_p_G + end end + """ $(TYPEDSIGNATURES) -Poisson bracket of a function and a Hamiltonian. +Internal implementation of Poisson bracket for autonomous, variable-dependent Hamiltonians. -Returns a Hamiltonian representing `{f, g}` where `g` is already a Hamiltonian. +Computes `{H, G}(x, p, v) = ∇ₚ H(x,p,v) · ∇ₓ G(x,p,v) - ∇ₓ H(x,p,v) · ∇ₚ G(x,p,v)`. -# Example -```julia-repl -julia> f = (x, p) -> x[2]^2 + 2x[1]^2 + p[1]^2 -julia> g = (x, p) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] -julia> G = Hamiltonian(g) -julia> Poisson(f, G)([1, 2], [2, 1]) # -20 - -julia> f = (t, x, p, v) -> t*v[1]*x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] -julia> g = (t, x, p, v) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] + t - v[2] -julia> G = Hamiltonian(g, autonomous=false, variable=true) -julia> Poisson(f, G)(2, [1, 2], [2, 1], [4, 4]) # -76 -``` +# Implementation + +Uses `ctgradient` to compute partial derivatives with respect to `x` and `p` separately. +The extra variable `v` is held constant during gradient computation. + +# Note + +This is an internal function. Users should call [`Poisson`](@ref) instead. + +See also: [`Poisson`](@ref), [`ctgradient`](@ref) """ -function Poisson( - f::Function, g::AbstractHamiltonian{TD,VD} -)::Hamiltonian where {TD<:TimeDependence,VD<:VariableDependence} - return Poisson(Hamiltonian(f, TD, VD), g) +function _Poisson(H::Function, G::Function, backend, ::Type{Autonomous}, ::Type{NonFixed}) + return function (x, p, v) + grad_x_H = ctgradient(y -> H(y, p, v), backend, x) + grad_p_H = ctgradient(q -> H(x, q, v), backend, p) + grad_x_G = ctgradient(y -> G(y, p, v), backend, x) + grad_p_G = ctgradient(q -> G(x, q, v), backend, p) + return grad_p_H' * grad_x_G - grad_x_H' * grad_p_G + end end + """ $(TYPEDSIGNATURES) -Poisson bracket of a Hamiltonian and a function. +Internal implementation of Poisson bracket for non-autonomous, time-fixed Hamiltonians. -Returns a Hamiltonian representing `{f, g}` where `f` is already a Hamiltonian. +Computes `{H, G}(t, x, p) = ∇ₚ H(t,x,p) · ∇ₓ G(t,x,p) - ∇ₓ H(t,x,p) · ∇ₚ G(t,x,p)`. -# Example -```julia-repl -julia> f = (x, p) -> x[2]^2 + 2x[1]^2 + p[1]^2 -julia> g = (x, p) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] -julia> F = Hamiltonian(f) -julia> Poisson(F, g)([1, 2], [2, 1]) # -20 - -julia> f = (t, x, p, v) -> t*v[1]*x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] -julia> g = (t, x, p, v) -> 3x[2]^2 - x[1]^2 + p[2]^2 + p[1] + t - v[2] -julia> F = Hamiltonian(f, autonomous=false, variable=true) -julia> Poisson(F, g)(2, [1, 2], [2, 1], [4, 4]) # -76 -``` +# Implementation + +Uses `ctgradient` to compute partial derivatives with respect to `x` and `p` separately. +Time `t` is held constant during gradient computation. + +# Note + +This is an internal function. Users should call [`Poisson`](@ref) instead. + +See also: [`Poisson`](@ref), [`ctgradient`](@ref) """ -function Poisson( - f::AbstractHamiltonian{TD,VD}, g::Function -)::Hamiltonian where {TD<:TimeDependence,VD<:VariableDependence} - return Poisson(f, Hamiltonian(g, TD, VD)) +function _Poisson(H::Function, G::Function, backend, ::Type{NonAutonomous}, ::Type{Fixed}) + return function (t, x, p) + grad_x_H = ctgradient(y -> H(t, y, p), backend, x) + grad_p_H = ctgradient(q -> H(t, x, q), backend, p) + grad_x_G = ctgradient(y -> G(t, y, p), backend, x) + grad_p_G = ctgradient(q -> G(t, x, q), backend, p) + return grad_p_H' * grad_x_G - grad_x_H' * grad_p_G + end end -# --------------------------------------------------------------------------- -# Macros -# @Lie [X, Y], for Lie brackets -# @Lie {X, Y}, for Poisson brackets """ $(TYPEDSIGNATURES) -Compute Lie or Poisson brackets. +Internal implementation of Poisson bracket for non-autonomous, variable-dependent Hamiltonians. -This macro provides a unified notation to define recursively nested Lie brackets (for vector fields) or Poisson brackets (for Hamiltonians). +Computes `{H, G}(t, x, p, v) = ∇ₚ H(t,x,p,v) · ∇ₓ G(t,x,p,v) - ∇ₓ H(t,x,p,v) · ∇ₚ G(t,x,p,v)`. -### Syntax +This is the most general case, handling time-dependent Hamiltonians with extra variables. -- `@Lie [F, G]`: computes the Lie bracket `[F, G]` of two vector fields. -- `@Lie [[F, G], H]`: supports arbitrarily nested Lie brackets. -- `@Lie {H, K}`: computes the Poisson bracket `{H, K}` of two Hamiltonians. -- `@Lie {{H, K}, L}`: supports arbitrarily nested Poisson brackets. -- `@Lie expr autonomous = false`: specifies a non-autonomous system. -- `@Lie expr variable = true`: indicates presence of an auxiliary variable `v`. +# Implementation -Keyword-like arguments can be provided to control the evaluation context for Poisson brackets with raw functions: -- `autonomous = Bool`: whether the system is time-independent (default: `true`). -- `variable = Bool`: whether the system depends on an extra variable `v` (default: `false`). +Uses `ctgradient` to compute partial derivatives with respect to `x` and `p` separately. +Time `t` and variable `v` are held constant during gradient computation. -### Bracket type detection +# Note -- Square brackets `[...]` denote Lie brackets between `VectorField` objects. -- Curly brackets `{...}` denote Poisson brackets between `Hamiltonian` objects or between raw functions. -- The macro automatically dispatches to `Lie` or `Poisson` depending on the input pattern. +This is an internal function. Users should call [`Poisson`](@ref) instead. -### Return +See also: [`Poisson`](@ref), [`ctgradient`](@ref) +""" +function _Poisson(H::Function, G::Function, backend, ::Type{NonAutonomous}, ::Type{NonFixed}) + return function (t, x, p, v) + grad_x_H = ctgradient(y -> H(t, y, p, v), backend, x) + grad_p_H = ctgradient(q -> H(t, x, q, v), backend, p) + grad_x_G = ctgradient(y -> G(t, y, p, v), backend, x) + grad_p_G = ctgradient(q -> G(t, x, q, v), backend, p) + return grad_p_H' * grad_x_G - grad_x_H' * grad_p_G + end +end -A callable object representing the specified Lie or Poisson bracket expression. The returned function can be evaluated like any other vector field or Hamiltonian. +# ============================================================================== +# Time Derivative +# ============================================================================== ---- +""" +$(TYPEDSIGNATURES) -### Examples +Partial derivative with respect to time of a function. -#### ■ Lie brackets with `VectorField` (autonomous) -```julia-repl -julia> F1 = VectorField(x -> [0, -x[3], x[2]]) -julia> F2 = VectorField(x -> [x[3], 0, -x[1]]) -julia> L = @Lie [F1, F2] -julia> L([1.0, 2.0, 3.0]) -3-element Vector{Float64}: - 2.0 - -1.0 - 0.0 -``` +Computes `∂f/∂t` for a function `f(t, ...)`. -#### ■ Lie brackets with `VectorField` (non-autonomous, with auxiliary variable) -```julia-repl -julia> F1 = VectorField((t, x, v) -> [0, -x[3], x[2]]; autonomous=false, variable=true) -julia> F2 = VectorField((t, x, v) -> [x[3], 0, -x[1]]; autonomous=false, variable=true) -julia> L = @Lie [F1, F2] -julia> L(0.0, [1.0, 2.0, 3.0], 1.0) -3-element Vector{Float64}: - 2.0 - -1.0 - 0.0 -``` +# Arguments -#### ■ Poisson brackets with `Hamiltonian` (autonomous) -```julia-repl -julia> H1 = Hamiltonian((x, p) -> x[1]^2 + p[2]^2) -julia> H2 = Hamiltonian((x, p) -> x[2]^2 + p[1]^2) -julia> P = @Lie {H1, H2} -julia> P([1.0, 1.0], [3.0, 2.0]) --4.0 -``` +- `f::Function`: A function with time as its first argument -#### ■ Poisson brackets with `Hamiltonian` (non-autonomous, with variable) -```julia-repl -julia> H1 = Hamiltonian((t, x, p, v) -> x[1]^2 + p[2]^2 + v; autonomous=false, variable=true) -julia> H2 = Hamiltonian((t, x, p, v) -> x[2]^2 + p[1]^2 + v; autonomous=false, variable=true) -julia> P = @Lie {H1, H2} -julia> P(1.0, [1.0, 3.0], [4.0, 2.0], 3.0) -8.0 -``` +# Returns -#### ■ Poisson brackets from raw functions -```julia-repl -julia> H1 = (x, p) -> x[1]^2 + p[2]^2 -julia> H2 = (x, p) -> x[2]^2 + p[1]^2 -julia> P = @Lie {H1, H2} -julia> P([1.0, 1.0], [3.0, 2.0]) --4.0 -``` +- A function `(t, args...) -> ∂f/∂t` that computes the partial derivative with respect to time + +# Example -#### ■ Poisson bracket with non-autonomous raw functions ```julia-repl -julia> H1 = (t, x, p) -> x[1]^2 + p[2]^2 + t -julia> H2 = (t, x, p) -> x[2]^2 + p[1]^2 + t -julia> P = @Lie {H1, H2} autonomous = false -julia> P(3.0, [1.0, 2.0], [4.0, 1.0]) --8.0 +julia> using CTFlows + +julia> f(t, x) = t^2 * x[1] + +julia> df = CTFlows.∂ₜ(f) + +julia> df(3.0, [2.0]) +12.0 ``` -#### ■ Nested brackets +See also: [`ctgradient`](@ref) +""" +∂ₜ(f) = (t, args...) -> ctgradient(y -> f(y, args...), t) + +# ============================================================================== +# @Lie Macro (V3: uses prefix system, no wrappers) +# ============================================================================== + +""" +$(TYPEDSIGNATURES) + +Macro for Lie brackets and Poisson brackets with mathematical notation. + +Uses type dispatch for better performance - expands to typed method calls. + +# Syntax + +- `@Lie [X, Y]` - Lie bracket (expands to `CTFlows.ad(X, Y, Autonomous, Fixed)`) +- `@Lie {H, G}` - Poisson bracket (expands to `CTFlows.Poisson(H, G, Autonomous, Fixed)`) +- `@Lie {H, G} autonomous=false` - With options (expands to `CTFlows.Poisson(H, G, NonAutonomous, Fixed)`) + +# Examples + ```julia-repl -julia> F = VectorField(x -> [-x[1], x[2], x[3]]) -julia> G = VectorField(x -> [x[3], -x[2], 0]) -julia> H = VectorField(x -> [0, 0, -x[1]]) -julia> nested = @Lie [[F, G], H] -julia> nested([1.0, 2.0, 3.0]) -3-element Vector{Float64}: +julia> using CTFlows + +julia> X(x) = [x[2], -x[1]] + +julia> Y(x) = [x[1], x[2]] + +julia> Z = CTFlows.@Lie [X, Y] + +julia> Z([1.0, 2.0]) +2-element Vector{Float64}: + -1.0 2.0 - 0.0 - -6.0 -``` -```julia-repl -julia> H1 = (x, p) -> x[2]*x[1]^2 + p[1]^2 -julia> H2 = (x, p) -> x[1]*p[2]^2 -julia> H3 = (x, p) -> x[1]*p[2] + x[2]*p[1] -julia> nested_poisson = @Lie {{H1, H2}, H3} -julia> nested_poisson([1.0, 2.0], [0.5, 1.0]) -14.0 -``` +julia> H(x, p) = x[1]^2 + p[1]^2 -#### ■ Mixed expressions with arithmetic -```julia-repl -julia> F1 = VectorField(x -> [0, -x[3], x[2]]) -julia> F2 = VectorField(x -> [x[3], 0, -x[1]]) -julia> x = [1.0, 2.0, 3.0] -julia> @Lie [F1, F2](x) + 3 * [F1, F2](x) -3-element Vector{Float64}: - 8.0 - -4.0 - 0.0 -``` +julia> G(x, p) = x[2]^2 + p[2]^2 -```julia-repl -julia> H1 = (x, p) -> x[1]^2 -julia> H2 = (x, p) -> p[1]^2 -julia> H3 = (x, p) -> x[1]*p[1] -julia> x = [1.0, 2.0, 3.0] -julia> p = [3.0, 2.0, 1.0] -julia> @Lie {H1, H2}(x, p) + 2 * {H2, H3}(x, p) -24.0 +julia> PB = CTFlows.@Lie {H, G} + +julia> PB([1.0, 2.0], [0.5, 0.5]) +0.0 ``` + +See also: [`ad`](@ref), [`Poisson`](@ref) """ macro Lie(expr::Expr, args...) - autonomous = true - variable = false + # Parse options - default values + autonomous = __autonomous() + variable = __variable() - # Parse keyword args for arg in args - @match arg begin - :(autonomous = $a) => (autonomous = a) - :(variable = $a) => (variable = a) - _ => throw(ArgumentError("Invalid argument: $arg")) + if @capture(arg, autonomous = val_) + autonomous = val + elseif @capture(arg, variable = val_) + variable = val end end - # Check for mixed usage of Lie and Poisson brackets - has_lie = Ref(false) - has_poisson = Ref(false) - - function check_mixed_usage(x) - function walker(e) - if @capture(e, [_, _]) - has_lie[] = true - elseif @capture(e, {_, _}) - has_poisson[] = true - end - return e - end - postwalk(walker, x) - - if has_lie[] && has_poisson[] - throw( - ArgumentError("Cannot mix Lie and Poisson brackets in the same expression.") - ) - end - return nothing - end + # Convert Bool to Type at compile time + TD = autonomous ? :Autonomous : :NonAutonomous + VD = variable ? :NonFixed : :Fixed - check_mixed_usage(expr) + prefix = diffgeo_prefix() - # Transform Lie and Poisson bracket expressions function fun(x) is_lie, is_poisson = @capture(x, [a_, b_]), @capture(x, {c_, d_}) if is_lie - # Just return Lie call with interpolation - return :(Lie($a, $b)) + # Lie bracket: @Lie [X, Y] -> CTFlows.ad(X, Y, CTFlows.Autonomous, CTFlows.Fixed) + return :( + $prefix.ad($a, $b, $prefix.$TD, $prefix.$VD) + ) elseif is_poisson - # Return a quoted block with if...else for runtime type checks - return quote - if isa($c, Function) && isa($d, Function) - Poisson( - Hamiltonian($c; autonomous=($(autonomous)), variable=($(variable))), - Hamiltonian($d; autonomous=($(autonomous)), variable=($(variable))), - ) - else - Poisson($c, $d) - end - end + # Poisson bracket: @Lie {H, G} -> CTFlows.Poisson(H, G, CTFlows.Autonomous, CTFlows.Fixed) + return :( + $prefix.Poisson($c, $d, $prefix.$TD, $prefix.$VD) + ) else return x end diff --git a/src/types.jl b/src/types.jl index 6a55bda..798c11a 100644 --- a/src/types.jl +++ b/src/types.jl @@ -103,118 +103,6 @@ struct Mayer{TF<:Function,VD<:VariableDependence} f::TF end -""" -$(TYPEDEF) -""" -abstract type AbstractHamiltonian{TD<:TimeDependence,VD<:VariableDependence} end - -""" -$(TYPEDEF) -""" -abstract type AbstractVectorField{TD<:TimeDependence,VD<:VariableDependence} end - -""" -$(TYPEDEF) - -Encodes the Hamiltonian function `H = ⟨p, f⟩ + L` in optimal control. - -# Fields -- `f`: a callable of the form: - - `f(x, p)` - - `f(t, x, p)` - - `f(x, p, v)` - - `f(t, x, p, v)` - -# Type Parameters -- `TD`: `Autonomous` or `NonAutonomous` -- `VD`: `Fixed` or `NonFixed` - -# Example -```julia-repl -julia> Hf(x, p) = dot(p, [x[2], -x[1]]) -julia> H = Hamiltonian{typeof(Hf), Autonomous, Fixed}(Hf) -julia> H([1.0, 0.0], [1.0, 1.0]) -``` -""" -struct Hamiltonian{TF<:Function,TD<:TimeDependence,VD<:VariableDependence} <: - AbstractHamiltonian{TD,VD} - f::TF -end - -""" -$(TYPEDEF) - -Represents a dynamical system `dx/dt = f(...)` as a vector field. - -# Fields -- `f`: a callable of the form: - - `f(x)` - - `f(t, x)` - - `f(x, v)` - - `f(t, x, v)` - -# Example -```julia-repl -f(x) = [x[2], -x[1]] -vf = VectorField{typeof(f), Autonomous, Fixed}(f) -vf([1.0, 0.0]) -``` -""" -struct VectorField{TF<:Function,TD<:TimeDependence,VD<:VariableDependence} <: - AbstractVectorField{TD,VD} - f::TF -end - -""" -$(TYPEDEF) - -Represents the Hamiltonian vector field associated to a Hamiltonian function, -typically defined as `(∂H/∂p, -∂H/∂x)`. - -# Fields -- `f`: a callable implementing the Hamiltonian vector field. - -# Example -```julia-repl -julia> f(x, p) = [p[2], -p[1], -x[1], -x[2]] -julia> XH = HamiltonianVectorField{typeof(f), Autonomous, Fixed}(f) -julia> XH([1.0, 0.0], [0.5, 0.5]) -``` -""" -struct HamiltonianVectorField{TF<:Function,TD<:TimeDependence,VD<:VariableDependence} <: - AbstractVectorField{TD,VD} - f::TF -end - -""" -$(TYPEDEF) - -Lifts a vector field `X` into a Hamiltonian function using the canonical symplectic structure. - -This is useful to convert a vector field into a Hamiltonian via the identity: -`H(x, p) = ⟨p, X(x)⟩`. - -# Constructor -Use `HamiltonianLift(X::VectorField)` where `X` is a `VectorField{...}`. - -# Example -```julia-repl -f(x) = [x[2], -x[1]] -julia> X = VectorField{typeof(f), Autonomous, Fixed}(f) -julia> H = HamiltonianLift(X) -julia> H([1.0, 0.0], [0.5, 0.5]) -``` -""" -struct HamiltonianLift{TV<:VectorField,TD<:TimeDependence,VD<:VariableDependence} <: - AbstractHamiltonian{TD,VD} - X::TV - function HamiltonianLift( - X::VectorField{<:Function,TD,VD} - ) where {TD<:TimeDependence,VD<:VariableDependence} - return new{typeof(X),TD,VD}(X) - end -end - """ $(TYPEDEF) @@ -384,315 +272,6 @@ function (F::Mayer{<:Function,NonFixed})(x0::State, xf::State, v::Variable)::ctN end # -------------------------------------------------------------------------------------------------- -""" -$(TYPEDSIGNATURES) - -Construct a Hamiltonian function wrapper. - -- `f`: a function representing the Hamiltonian. -- `autonomous`: whether `f` is autonomous (default via `__autonomous()`). -- `variable`: whether `f` depends on an extra variable argument (default via `__variable()`). - -Returns a `Hamiltonian{TF, TD, VD}` callable struct. -""" -function Hamiltonian( - f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -) - TD = autonomous ? Autonomous : NonAutonomous - VD = variable ? NonFixed : Fixed - return Hamiltonian{typeof(f),TD,VD}(f) -end - -""" -$(TYPEDSIGNATURES) - -Construct a Hamiltonian function wrapper with explicit time and variable dependence types. -""" -function Hamiltonian( - f::Function, TD::Type{<:TimeDependence}, VD::Type{<:VariableDependence} -) - return Hamiltonian{typeof(f),TD,VD}(f) -end - -function (F::Hamiltonian{<:Function,Autonomous,Fixed})(x::State, p::Costate)::ctNumber - return F.f(x, p) -end - -function (F::Hamiltonian{<:Function,Autonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return F.f(x, p) -end - -function (F::Hamiltonian{<:Function,Autonomous,NonFixed})( - x::State, p::Costate, v::Variable -)::ctNumber - return F.f(x, p, v) -end - -function (F::Hamiltonian{<:Function,Autonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return F.f(x, p, v) -end - -function (F::Hamiltonian{<:Function,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate -)::ctNumber - return F.f(t, x, p) -end - -function (F::Hamiltonian{<:Function,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return F.f(t, x, p) -end - -function (F::Hamiltonian{<:Function,NonAutonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return F.f(t, x, p, v) -end - -# --------------------------------------------------------------------------- -""" -$(TYPEDSIGNATURES) - -Construct a HamiltonianLift from a vector field function. - -- `f`: function defining the vector field. -- `autonomous`: whether `f` is autonomous. -- `variable`: whether `f` depends on an extra variable argument. - -Returns a `HamiltonianLift` wrapping the vector field. -""" -function HamiltonianLift( - f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -) - TD = autonomous ? Autonomous : NonAutonomous - VD = variable ? NonFixed : Fixed - return HamiltonianLift(VectorField(f, TD, VD)) -end - -""" -$(TYPEDSIGNATURES) - -Construct a HamiltonianLift with explicit time and variable dependence types. -""" -function HamiltonianLift( - f::Function, TD::Type{<:TimeDependence}, VD::Type{<:VariableDependence} -) - return HamiltonianLift(VectorField(f, TD, VD)) -end - -function (H::HamiltonianLift{<:VectorField,Autonomous,Fixed})( - x::State, p::Costate -)::ctNumber - return p' * H.X(x) -end - -function (H::HamiltonianLift{<:VectorField,Autonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return p' * H.X(x) -end - -function (H::HamiltonianLift{<:VectorField,Autonomous,NonFixed})( - x::State, p::Costate, v::Variable -)::ctNumber - return p' * H.X(x, v) -end - -function (H::HamiltonianLift{<:VectorField,Autonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return p' * H.X(x, v) -end - -function (H::HamiltonianLift{<:VectorField,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate -)::ctNumber - return p' * H.X(t, x) -end - -function (H::HamiltonianLift{<:VectorField,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return p' * H.X(t, x) -end - -function (H::HamiltonianLift{<:VectorField,NonAutonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::ctNumber - return p' * H.X(t, x, v) -end - -# -------------------------------------------------------------------------------------------------- -""" -$(TYPEDSIGNATURES) - -Construct a Hamiltonian vector field from a function `f`. - -- `autonomous`: whether `f` is autonomous. -- `variable`: whether `f` depends on an extra variable argument. - -Returns a `HamiltonianVectorField{TF, TD, VD}` callable struct. -""" -function HamiltonianVectorField( - f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -) - TD = autonomous ? Autonomous : NonAutonomous - VD = variable ? NonFixed : Fixed - return HamiltonianVectorField{typeof(f),TD,VD}(f) -end - -""" -$(TYPEDSIGNATURES) - -Construct a Hamiltonian vector field with explicit time and variable dependence. -""" -function HamiltonianVectorField( - f::Function, TD::Type{<:TimeDependence}, VD::Type{<:VariableDependence} -) - return HamiltonianVectorField{typeof(f),TD,VD}(f) -end - -function (F::HamiltonianVectorField{<:Function,Autonomous,Fixed})( - x::State, p::Costate -)::Tuple{DState,DCostate} - return F.f(x, p) -end - -function (F::HamiltonianVectorField{<:Function,Autonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::Tuple{DState,DCostate} - return F.f(x, p) -end - -function (F::HamiltonianVectorField{<:Function,Autonomous,NonFixed})( - x::State, p::Costate, v::Variable -)::Tuple{DState,DCostate} - return F.f(x, p, v) -end - -function (F::HamiltonianVectorField{<:Function,Autonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::Tuple{DState,DCostate} - return F.f(x, p, v) -end - -function (F::HamiltonianVectorField{<:Function,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate -)::Tuple{DState,DCostate} - return F.f(t, x, p) -end - -function (F::HamiltonianVectorField{<:Function,NonAutonomous,Fixed})( - t::Time, x::State, p::Costate, v::Variable -)::Tuple{DState,DCostate} - return F.f(t, x, p) -end - -function (F::HamiltonianVectorField{<:Function,NonAutonomous,NonFixed})( - t::Time, x::State, p::Costate, v::Variable -)::Tuple{DState,DCostate} - return F.f(t, x, p, v) -end - -# -------------------------------------------------------------------------------------------------- -""" -$(TYPEDSIGNATURES) - -Create a `VectorField` object wrapping the function `f`. - -# Arguments -- `f::Function`: The vector field function. -- `autonomous::Bool` (optional): If `true`, the vector field is autonomous (time-independent). Defaults to `__autonomous()`. -- `variable::Bool` (optional): If `true`, the vector field depends on control or decision variables (non-fixed). Defaults to `__variable()`. - -# Returns -- A `VectorField{typeof(f),TD,VD}` object where `TD` encodes time dependence and `VD` encodes variable dependence. - -# Details -The `VectorField` object can be called with different argument signatures depending on the time and variable dependence. - -# Examples -```julia-repl -julia> f(x) = [-x[2], x[1]] -julia> vf = VectorField(f, autonomous=true, variable=false) -julia> vf([1.0, 0.0]) # returns [-0.0, 1.0] -``` -""" -function VectorField( - f::Function; autonomous::Bool=__autonomous(), variable::Bool=__variable() -) - TD = autonomous ? Autonomous : NonAutonomous - VD = variable ? NonFixed : Fixed - return VectorField{typeof(f),TD,VD}(f) -end - -""" -$(TYPEDSIGNATURES) - -Create a `VectorField` object with explicit time and variable dependence types. - -# Arguments -- `f::Function`: The vector field function. -- `TD`: Type indicating time dependence (`Autonomous` or `NonAutonomous`). -- `VD`: Type indicating variable dependence (`Fixed` or `NonFixed`). - -# Returns -- A `VectorField{typeof(f),TD,VD}` object. -""" -function VectorField( - f::Function, TD::Type{<:TimeDependence}, VD::Type{<:VariableDependence} -) - return VectorField{typeof(f),TD,VD}(f) -end - -# Dispatch for autonomous fixed variable dependence: f(x) -function (F::VectorField{<:Function,Autonomous,Fixed})(x::State)::ctVector - return F.f(x) -end - -# Dispatch for autonomous fixed variable dependence with unused args: f(x) -function (F::VectorField{<:Function,Autonomous,Fixed})( - t::Time, x::State, v::Variable -)::ctVector - return F.f(x) -end - -# Dispatch for autonomous non-fixed variable dependence: f(x, v) -function (F::VectorField{<:Function,Autonomous,NonFixed})(x::State, v::Variable)::ctVector - return F.f(x, v) -end - -# Dispatch for autonomous non-fixed variable dependence with time arg: f(x, v) -function (F::VectorField{<:Function,Autonomous,NonFixed})( - t::Time, x::State, v::Variable -)::ctVector - return F.f(x, v) -end - -# Dispatch for non-autonomous fixed variable dependence: f(t, x) -function (F::VectorField{<:Function,NonAutonomous,Fixed})(t::Time, x::State)::ctVector - return F.f(t, x) -end - -# Dispatch for non-autonomous fixed variable dependence with unused variable: f(t, x) -function (F::VectorField{<:Function,NonAutonomous,Fixed})( - t::Time, x::State, v::Variable -)::ctVector - return F.f(t, x) -end - -# Dispatch for non-autonomous non-fixed variable dependence: f(t, x, v) -function (F::VectorField{<:Function,NonAutonomous,NonFixed})( - t::Time, x::State, v::Variable -)::ctVector - return F.f(t, x, v) -end - # -------------------------------------------------------------------------------------------------- """ $(TYPEDSIGNATURES) diff --git a/src/utils.jl b/src/utils.jl index f3f958d..8eb33bc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,6 +5,7 @@ Compute the derivative of a scalar function `f` at a scalar point `x`. # Arguments - `f::Function`: A scalar-valued function. +- `backend`: Automatic differentiation backend (default: `__backend()` = `AutoForwardDiff()`) - `x::ctNumber`: A scalar input. # Returns @@ -15,109 +16,60 @@ Compute the derivative of a scalar function `f` at a scalar point `x`. julia> ctgradient(x -> x^2, 3.0) # returns 6.0 ``` """ -function ctgradient(f::Function, x::ctNumber) - return ForwardDiff.derivative(x -> f(x), x) -end - -""" -$(TYPEDSIGNATURES) - -Compute the gradient of a scalar function `f` at a vector point `x`. - -# Arguments -- `f::Function`: A scalar-valued function accepting a vector input. -- `x`: A vector of numbers. - -# Returns -- A vector representing the gradient ∇f(x). - -# Example -```julia-repl -julia> ctgradient(x -> sum(x.^2), [1.0, 2.0]) # returns [2.0, 4.0] -``` -""" -function ctgradient(f::Function, x) - return ForwardDiff.gradient(f, x) +function ctgradient(f::Function, backend::AbstractADType, x::ctNumber) + return derivative(f, backend, x) end """ $(TYPEDSIGNATURES) -Compute the gradient of a `VectorField` at a given point. - -# Arguments -- `X::VectorField`: A vector field object with a callable function `X.f`. -- `x`: A scalar or vector input. - -# Returns -- The derivative or gradient depending on the type of `x`. - -# Example -```julia-repl -julia> X = VectorField(x -> x^2) -julia> ctgradient(X, 2.0) # returns 4.0 -``` -""" -ctgradient(X::VectorField, x) = ctgradient(X.f, x) - -""" -$(TYPEDSIGNATURES) - -Compute the Jacobian of a vector-valued function `f` at a scalar point `x`. +Compute the derivative of a scalar function `f` at a scalar point `x` using default backend. # Arguments -- `f::Function`: A vector-valued function. +- `f::Function`: A scalar-valued function. - `x::ctNumber`: A scalar input. # Returns -- A matrix representing the Jacobian Jf(x). - -# Example -```julia-repl -julia> f(x) = [sin(x), cos(x)] -julia> ctjacobian(f, 0.0) # returns a 2×1 matrix -``` +- The derivative of `f` evaluated at `x`. """ -function ctjacobian(f::Function, x::ctNumber) - return ForwardDiff.jacobian(x -> f(x[1]), [x]) +function ctgradient(f::Function, x::ctNumber) + return ctgradient(f, __backend(), x) end """ $(TYPEDSIGNATURES) -Compute the Jacobian of a vector-valued function `f` at a vector point `x`. +Compute the gradient of a scalar function `f` at a vector point `x`. # Arguments -- `f::Function`: A vector-valued function. -- `x`: A vector input. +- `f::Function`: A scalar-valued function accepting a vector input. +- `backend`: Automatic differentiation backend (default: `__backend()` = `AutoForwardDiff()`) +- `x`: A vector of numbers. # Returns -- A matrix representing the Jacobian Jf(x). +- A vector representing the gradient ∇f(x). # Example ```julia-repl -julia> f(x) = [x[1]^2, x[2]^2] -julia> ctjacobian(f, [1.0, 2.0]) # returns [2.0 0.0; 0.0 4.0] +julia> ctgradient(x -> sum(x.^2), [1.0, 2.0]) # returns [2.0, 4.0] ``` """ -ctjacobian(f::Function, x) = ForwardDiff.jacobian(f, x) +function ctgradient(f::Function, backend::AbstractADType, x) + return gradient(f, backend, x) +end """ $(TYPEDSIGNATURES) -Compute the Jacobian of a `VectorField` at a given point. +Compute the gradient of a scalar function `f` at a vector point `x` using default backend. # Arguments -- `X::VectorField`: A vector field object with a callable function `X.f`. -- `x`: A scalar or vector input. +- `f::Function`: A scalar-valued function accepting a vector input. +- `x`: A vector of numbers. # Returns -- A matrix representing the Jacobian of `X` at `x`. - -# Example -```julia-repl -julia> X = VectorField(x -> [x[1]^2, x[2]]) -julia> ctjacobian(X, [1.0, 3.0]) # returns [2.0 0.0; 0.0 1.0] -``` +- A vector representing the gradient ∇f(x). """ -ctjacobian(X::VectorField, x) = ctjacobian(X.f, x) +function ctgradient(f::Function, x) + return ctgradient(f, __backend(), x) +end diff --git a/test/test_default.jl b/test/test_default.jl index bbb4b00..16ee25c 100644 --- a/test/test_default.jl +++ b/test/test_default.jl @@ -1,11 +1,29 @@ function test_default() - Test.@test CTFlowsODE.__abstol() isa Real - Test.@test CTFlowsODE.__abstol() > 0 - Test.@test CTFlowsODE.__abstol() < 1 - Test.@test CTFlowsODE.__reltol() isa Real - Test.@test CTFlowsODE.__reltol() > 0 - Test.@test CTFlowsODE.__reltol() < 1 - Test.@test CTFlowsODE.__saveat() isa Vector - Test.@test CTFlowsODE.__alg() isa Tsit5 - Test.@test CTFlowsODE.__tstops() isa Vector{<:CTFlows.Time} + # Test core default functions + Test.@testset "Core defaults" begin + # Test __autonomous() + Test.@test CTFlows.__autonomous() isa Bool + Test.@test CTFlows.__autonomous() == true + + # Test __variable() + Test.@test CTFlows.__variable() isa Bool + Test.@test CTFlows.__variable() == false + + # Test __backend() + backend = CTFlows.__backend() + Test.@test backend isa DifferentiationInterface.AutoForwardDiff + end + + # Test extension defaults (CTFlowsODE) + Test.@testset "CTFlowsODE defaults" begin + Test.@test CTFlowsODE.__abstol() isa Real + Test.@test CTFlowsODE.__abstol() > 0 + Test.@test CTFlowsODE.__abstol() < 1 + Test.@test CTFlowsODE.__reltol() isa Real + Test.@test CTFlowsODE.__reltol() > 0 + Test.@test CTFlowsODE.__reltol() < 1 + Test.@test CTFlowsODE.__saveat() isa Vector + Test.@test CTFlowsODE.__alg() isa Tsit5 + Test.@test CTFlowsODE.__tstops() isa Vector{<:CTFlows.Time} + end end diff --git a/test/test_differential_geometry.jl b/test/test_differential_geometry.jl index 8733920..02a851d 100644 --- a/test/test_differential_geometry.jl +++ b/test/test_differential_geometry.jl @@ -1,818 +1,449 @@ -function test_differential_geometry() - - # - Hamiltonian = CTFlows.Hamiltonian - VectorField = CTFlows.VectorField - HamiltonianLift = CTFlows.HamiltonianLift - Lift = CTFlows.Lift - Lie = CTFlows.Lie - ⅋ = CTFlows.:(⅋) - ⋅ = CTFlows.:(⋅) - Poisson = CTFlows.Poisson - - # - ∅ = Vector{Real}() - dummy_function() = nothing - - @testset "Lifts" begin - @testset "HamiltonianLift from VectorField" begin - HL = HamiltonianLift( - VectorField(x -> [x[1]^2, x[2]^2]; autonomous=true, variable=false) - ) - @test HL([1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - VectorField((x, v) -> [x[1]^2, x[2]^2 + v]; autonomous=true, variable=true) - ) - @test HL([1, 0], [0, 1], 1) == 1 - @test HL(1, [1, 0], [0, 1], 1) == 1 - HL = HamiltonianLift( - VectorField( - (t, x) -> [t + x[1]^2, x[2]^2]; autonomous=false, variable=false - ), - ) - @test HL(1, [1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - VectorField( - (t, x, v) -> [t + x[1]^2, x[2]^2 + v]; autonomous=false, variable=true - ), - ) - @test HL(1, [1, 0], [0, 1], 1) == 1 - end - - @testset "HamiltonianLift from Function" begin - HL = HamiltonianLift(x -> [x[1]^2, x[2]^2]; autonomous=true, variable=false) - @test HL([1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - (x, v) -> [x[1]^2, x[2]^2 + v]; autonomous=true, variable=true - ) - @test HL([1, 0], [0, 1], 1) == 1 - @test HL(1, [1, 0], [0, 1], 1) == 1 - HL = HamiltonianLift( - (t, x) -> [t + x[1]^2, x[2]^2]; autonomous=false, variable=false - ) - @test HL(1, [1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - (t, x, v) -> [t + x[1]^2, x[2]^2 + v]; autonomous=false, variable=true - ) - @test HL(1, [1, 0], [0, 1], 1) == 1 - HL = HamiltonianLift(x -> [x[1]^2, x[2]^2]; autonomous=true, variable=false) - @test HL([1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - (x, v) -> [x[1]^2, x[2]^2 + v]; autonomous=true, variable=true - ) - @test HL([1, 0], [0, 1], 1) == 1 - @test HL(1, [1, 0], [0, 1], 1) == 1 - HL = HamiltonianLift( - (t, x) -> [t + x[1]^2, x[2]^2]; autonomous=false, variable=false - ) - @test HL(1, [1, 0], [0, 1]) == 0 - @test HL(1, [1, 0], [0, 1], 1) == 0 - HL = HamiltonianLift( - (t, x, v) -> [t + x[1]^2, x[2]^2 + v]; autonomous=false, variable=true - ) - @test HL(1, [1, 0], [0, 1], 1) == 1 - end - - @testset "from VectorFields" begin - - # autonomous case - X = VectorField(x -> 2x) - H = Lift(X) - Test.@test H(1, 1) == 2 - Test.@test H([1, 2], [3, 4]) == 22 - Test.@test_throws MethodError H([1, 2], 1) - Test.@test_throws MethodError H(1, [1, 2]) - - # nonautonomous case - X = VectorField((t, x) -> 2x; autonomous=false, variable=false) - H = Lift(X) - Test.@test H(1, 1, 1) == 2 - Test.@test H(1, [1, 2], [3, 4]) == 22 - - # autonomous nonfixed case - X = VectorField((x, v) -> 2x; autonomous=true, variable=true) - H = Lift(X) - Test.@test H(1, 1, 1) == 2 - Test.@test H([1, 2], [3, 4], 1) == 22 - Test.@test_throws MethodError H([1, 2], 1) - Test.@test_throws MethodError H(1, [1, 2]) - - # nonautonomous nonfixed case - X = VectorField((t, x, v) -> 2x; autonomous=false, variable=true) - H = Lift(X) - Test.@test H(1, 1, 1, 1) == 2 - Test.@test H(1, [1, 2], [3, 4], 1) == 22 - end - - @testset "from Function" begin - - # autonomous case - X::Function = x -> 2x - H = Lift(X) - Test.@test H(1, 1) == 2 - Test.@test H([1, 2], [3, 4]) == 22 - - # nonautonomous case - Xt::Function = (t, x) -> 2x - H = Lift(Xt; autonomous=false) - Test.@test H(1, 1, 1) == 2 - Test.@test H(1, [1, 2], [3, 4]) == 22 - - # autonomous nonfixed case - Xv::Function = (x, v) -> 2x - H = Lift(Xv; variable=true) - Test.@test H(1, 1, 1) == 2 - Test.@test H([1, 2], [3, 4], 1) == 22 - - # nonautonomous nonfixed case - Xtv::Function = (t, x, v) -> 2x - H = Lift(Xtv; autonomous=false, variable=true) - Test.@test H(1, 1, 1, 1) == 2 - Test.@test H(1, [1, 2], [3, 4], 1) == 22 - - # overload - F::Function = x -> 2x - H_F(x, p) = Lift(F)(x, p) - H_F(y) = H_F(y, y) - Test.@test H_F(1) == 2 - end - end # tests for Lift - - @testset "Directional derivative of a scalar function" begin - - # autonomous, dim 2 - φ = x -> [x[2], -x[1]] - X = VectorField(φ) - f = x -> x[1]^2 + x[2]^2 - Test.@test (X ⋅ f)([1, 2]) == 0 - Test.@test (X ⋅ f)([1, 2]) == Lie(X, f)([1, 2]) - Test.@test Lie(φ, f)([1, 2]) == 0 - - # autonomous, dim 1 - φ = x -> 2x - X = VectorField(φ) - f = x -> x^2 - Test.@test (X ⋅ f)(1) == 4 - Test.@test (φ ⋅ f)(1) == 4 - Test.@test Lie(φ, f)(1) == 4 - - # nonautonomous, dim 2 - φ = (t, x) -> [t + x[2], -x[1]] - X = VectorField(φ; autonomous=false, variable=false) - f = (t, x) -> t + x[1]^2 + x[2]^2 - Test.@test (X ⋅ f)(1, [1, 2]) == 2 - Test.@test Lie(φ, f; autonomous=false, variable=false)(1, [1, 2]) == 2 - - # nonautonomous, dim 1 - φ = (t, x) -> 2x + t - X = VectorField(φ; autonomous=false, variable=false) - f = (t, x) -> t + x^2 - Test.@test (X ⋅ f)(1, 1) == 6 - Test.@test Lie(φ, f; autonomous=false, variable=false)(1, 1) == 6 - - # autonomous, nonfixed, dim 2 - φ = (x, v) -> [x[2] + v[1], -x[1] + v[2]] - X = VectorField(φ; autonomous=true, variable=true) - f = (x, v) -> x[1]^2 + x[2]^2 - Test.@test (X ⋅ f)([1, 2], [2, 1]) == 8 - Test.@test Lie(φ, f; autonomous=true, variable=true)([1, 2], [2, 1]) == 8 - - # autonomous, nonfixed, dim 1 - φ = (x, v) -> 2x + v - X = VectorField(φ; autonomous=true, variable=true) - f = (x, v) -> x^2 - Test.@test (X ⋅ f)(1, 1) == 6 - Test.@test Lie(φ, f; autonomous=true, variable=true)(1, 1) == 6 - - # nonautonomous, nonfixed, dim 2 - φ = (t, x, v) -> [t + x[2] + v[1], -x[1] + v[2]] - X = VectorField(φ; autonomous=false, variable=true) - f = (t, x, v) -> t + x[1]^2 + x[2]^2 - Test.@test (X ⋅ f)(1, [1, 2], [2, 1]) == 10 - Test.@test Lie(φ, f; autonomous=false, variable=true)(1, [1, 2], [2, 1]) == 10 - - # nonautonomous, nonfixed, dim 1 - φ = (t, x, v) -> 2x + t + v - X = VectorField(φ; autonomous=false, variable=true) - f = (t, x, v) -> t + x^2 - Test.@test (X ⋅ f)(1, 1, 1) == 8 - Test.@test Lie(φ, f; autonomous=false, variable=true)(1, 1, 1) == 8 +using CTFlows +using Test +using LinearAlgebra + +# Helper for testing +function isapprox_vec(a, b; atol=1e-6) + return norm(a - b) < atol +end + +@testset "Differential Geometry V3" begin + + @testset "ad() - Lie Derivative" begin + # Autonomous + X(x) = [x[2], -x[1]] + f(x) = x[1]^2 + x[2]^2 + Lf = CTFlows.ad(X, f) + @test Lf([1.0, 2.0]) ≈ 0.0 atol = 1e-6 + + # Non-autonomous + X_na(t, x) = [t * x[2], -x[1]] + f_na(t, x) = t + x[1]^2 + Lf_na = CTFlows.ad(X_na, f_na; autonomous=false) + # Lf(t, x) = ∂f/∂x * X + ∂f/∂t (Wait, Lie derivative usually only space derivative?) + # Standard definition L_X f = ∇f · X. Time derivative is separate usually. + # Implementation check: + # ad uses derivative(g, ...). g(s) = foo(t, x + s*X) + # d/ds g(0) = ∇x foo · X. Correct. + + # ∂f/∂x = [2x1, 0] + # X = [t*x2, -x1] + # dot = 2x1*t*x2 + # t=2, x=[1, 2] -> 2*1*2*2 = 8 + @test Lf_na(2.0, [1.0, 2.0]) ≈ 8.0 atol = 1e-6 end - @testset "Directional derivative of a vector field" begin - - # autonomous, dim 2 - X = VectorField(x -> [x[2], -x[1]]) - Y = VectorField(x -> [x[1], x[2]]) - Test.@test (X ⅋ Y)([1, 2]) == [2, -1] - - # autonomous, dim 1 - X = VectorField(x -> 2x) - Y = VectorField(x -> 3x) - Test.@test (X ⅋ Y)(1) == 6 - - # nonautonomous, dim 2 - X = VectorField((t, x) -> [t + x[2], -x[1]]; autonomous=false, variable=false) - Y = VectorField((t, x) -> [t + x[1], x[2]]; autonomous=false, variable=false) - Test.@test (X ⅋ Y)(1, [1, 2]) == [3, -1] - - # nonautonomous, dim 1 - X = VectorField((t, x) -> 2x + t; autonomous=false, variable=false) - Y = VectorField((t, x) -> 3x + t; autonomous=false, variable=false) - Test.@test (X ⅋ Y)(1, 1) == 9 - - # autonomous, nonfixed, dim 1 - X = VectorField((x, v) -> 2x + v; autonomous=true, variable=true) - Y = VectorField((x, v) -> 3x + v; autonomous=true, variable=true) - Test.@test (X ⅋ Y)(1, 1) == 9 - - # nonautonomous, nonfixed, dim 1 - X = VectorField((t, x, v) -> t + 2x + v; autonomous=false, variable=true) - Y = VectorField((t, x, v) -> t + 3x + v; autonomous=false, variable=true) - Test.@test (X ⅋ Y)(1, 1, 1) == 12 - - # autonomous, nonfixed, dim 2 - X = VectorField( - (x, v) -> [v[1] + v[2] + x[2], -x[1]]; autonomous=true, variable=true - ) - Y = VectorField( - (x, v) -> [v[1] + v[2] + x[1], x[2]]; autonomous=true, variable=true - ) - Test.@test (X ⅋ Y)([1, 2], [2, 3]) == [7, -1] - - # nonautonomous, nonfixed, dim 2 - X = VectorField( - (t, x, v) -> [t + v[1] + v[2] + x[2], -x[1]]; autonomous=false, variable=true - ) - Y = VectorField( - (t, x, v) -> [v[1] + v[2] + x[1], x[2]]; autonomous=false, variable=true - ) - Test.@test (X ⅋ Y)(1, [1, 2], [2, 3]) == [8, -1] + @testset "ad() - Lie Bracket" begin + # Autonomous + X(x) = [x[2], 0.0] + Y(x) = [0.0, x[1]] + # [X, Y] = J_Y*X - J_X*Y + # J_Y = [0 0; 1 0], X = [x2, 0] -> J_Y*X = [0, x2] + # J_X = [0 1; 0 0], Y = [0, x1] -> J_X*Y = [x1, 0] + # [X, Y] = [0, x2] - [x1, 0] = [-x1, x2] + XY = CTFlows.ad(X, Y) + @test isapprox_vec(XY([1.0, 2.0]), [-1.0, 2.0]) + + # Macro check + XY_macro = CTFlows.@Lie [X, Y] + @test isapprox_vec(XY_macro([1.0, 2.0]), [-1.0, 2.0]) end - @testset "Lie bracket - length 2" begin - @testset "autonomous case" begin - f = x -> [x[2], 2x[1]] - g = x -> [3x[2], -x[1]] - X = VectorField(f) - Y = VectorField(g) - Test.@test Lie(X, Y)([1, 2]) == [7, -14] - end - - @testset "nonautonomous case" begin - f = (t, x) -> [t + x[2], -2x[1]] - g = (t, x) -> [t + 3x[2], -x[1]] - X = VectorField(f; autonomous=false, variable=false) - Y = VectorField(g; autonomous=false, variable=false) - Test.@test Lie(X, Y)(1, [1, 2]) == [-5, 11] - end - - @testset "autonomous nonfixed case" begin - f = (x, v) -> [x[2] + v, 2x[1]] - g = (x, v) -> [3x[2], v - x[1]] - X = VectorField(f; autonomous=true, variable=true) - Y = VectorField(g; variable=true) - Test.@test Lie(X, Y)([1, 2], 1) == [6, -15] - end - - @testset "nonautonomous nonfixed case" begin - f = (t, x, v) -> [t + x[2] + v, -2x[1] - v] - g = (t, x, v) -> [t + 3x[2] + v, -x[1] - v] - X = VectorField(f; autonomous=false, variable=true) - Y = VectorField(g; autonomous=false, variable=true) - Test.@test Lie(X, Y)(1, [1, 2], 1) == [-7, 12] - end - - @testset "mri example" begin - Γ = 2 - γ = 1 - δ = γ - Γ - F0 = VectorField(x -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]) - F1 = VectorField(x -> [0, -x[3], x[2]]) - F2 = VectorField(x -> [x[3], 0, -x[1]]) - F01 = Lie(F0, F1) - F02 = Lie(F0, F2) - F12 = Lie(F1, F2) - x = [1, 2, 3] - Test.@test F01(x) ≈ -[0, γ - δ * x[3], -δ * x[2]] atol = 1e-6 - Test.@test F02(x) ≈ -[-γ + δ * x[3], 0, δ * x[1]] atol = 1e-6 - Test.@test F12(x) ≈ -[-x[2], x[1], 0] atol = 1e-6 - end - - @testset "intrinsic definition" begin - X = VectorField(x -> [x[2]^2, -2x[1] * x[2]]) - Y = VectorField(x -> [x[1] * (1 + x[2]), 3x[2]^3]) - f = x -> x[1]^4 + 2x[2]^3 - x = [1, 2] - - # check real intrinsic order 2 definition of Lie bracket - Test.@test ((CTFlows.@Lie [X, Y]) ⋅ f)(x) == - ((X ⋅ (Y ⋅ f))(x) - (Y ⋅ (X ⋅ f))(x)) - end + @testset "Lift()" begin + f(x) = [x[1], x[2]] + H = CTFlows.Lift(f) + # H(x, p) = p' * f(x) = p1*x1 + p2*x2 + @test H([1.0, 2.0], [3.0, 4.0]) ≈ 1.0 * 3.0 + 2.0 * 4.0 atol = 1e-6 + + # Non-autonomous + g(t, x) = [t * x[1], x[2]] + Hg = CTFlows.Lift(g; autonomous=false) + @test Hg(2.0, [1.0, 2.0], [3.0, 4.0]) ≈ 3.0 * (2.0 * 1.0) + 4.0 * 2.0 atol = 1e-6 + end + + @testset "Poisson()" begin + H(x, p) = 0.5 * (p[1]^2 + p[2]^2) + G(x, p) = x[1] + # {H, G} = ∂H/∂p * ∂G/∂x - ∂H/∂x * ∂G/∂p + # ∂H/∂p = [p1, p2], ∂G/∂x = [1, 0] -> dot = p1 + # ∂H/∂x = [0, 0] -> term 2 is 0 + PB = CTFlows.Poisson(H, G) + @test PB([1.0, 2.0], [3.0, 4.0]) ≈ 3.0 atol = 1e-6 + + # Macro check + PB_macro = CTFlows.@Lie {H, G} + @test PB_macro([1.0, 2.0], [3.0, 4.0]) ≈ 3.0 atol = 1e-6 + + # Scalar case + H_scalar(x, p) = 0.5 * (p^2 + x^2) + G_scalar(x, p) = x + # {H, G} = ∂H/∂p * ∂G/∂x - ∂H/∂x * ∂G/∂p + # ∂H/∂p = p, ∂G/∂x = 1 -> p*1 = p + # ∂H/∂x = x, ∂G/∂p = 0 -> x*0 = 0 + # {H, G} = p + PB_scalar = CTFlows.Poisson(H_scalar, G_scalar) + @test PB_scalar(1.0, 3.0) ≈ 3.0 atol = 1e-6 end - @testset "Poisson bracket" begin - @testset "autonomous case" begin - f = (x, p) -> x[2]^2 + 2x[1]^2 + p[1]^2 - g = (x, p) -> 3x[2]^2 + -x[1]^2 + p[2]^2 + p[1] - h = (x, p) -> x[2]^2 + -2x[1]^2 + p[1]^2 - 2p[2]^2 - f₊g = (x, p) -> f(x, p) + g(x, p) - fg = (x, p) -> f(x, p) * g(x, p) - # - F = Hamiltonian(f) - G = Hamiltonian(g) - H = Hamiltonian(h) - F₊G = Hamiltonian(f₊g) - FG = Hamiltonian(fg) - # - Test.@test Poisson(f, g)([1, 2], [2, 1]) == -20 - Test.@test Poisson(f, G)([1, 2], [2, 1]) == -20 - Test.@test Poisson(F, g)([1, 2], [2, 1]) == -20 - # - Test.@test Poisson(F, Hamiltonian((x, p) -> 42))([1, 2], [2, 1]) == 0 - Test.@test Poisson(F, G)([1, 2], [2, 1]) == -20 - Test.@test Poisson(F, G)([1, 2], [2, 1]) == -Poisson(G, F)([1, 2], [2, 1]) # anticommutativity - Test.@test Poisson(F₊G, H)([1, 2], [2, 1]) == - Poisson(F, H)([1, 2], [2, 1]) + Poisson(G, H)([1, 2], [2, 1]) # bilinearity 1 - Test.@test Poisson(H, F₊G)([1, 2], [2, 1]) == - Poisson(H, F)([1, 2], [2, 1]) + Poisson(H, G)([1, 2], [2, 1]) # bilinearity 2 - Test.@test Poisson(FG, H)([1, 2], [2, 1]) == - Poisson(F, H)([1, 2], [2, 1]) * G([1, 2], [2, 1]) + - F([1, 2], [2, 1]) * Poisson(G, H)([1, 2], [2, 1]) # Liebniz's rule - Test.@test Poisson(F, Poisson(G, H))([1, 2], [2, 1]) + - Poisson(G, Poisson(H, F))([1, 2], [2, 1]) + - Poisson(H, Poisson(F, G))([1, 2], [2, 1]) == 0 # Jacobi identity - end - - @testset "nonautonomous case" begin - f = (t, x, p) -> t * x[2]^2 + 2x[1]^2 + p[1]^2 + t - g = (t, x, p) -> 3x[2]^2 + -x[1]^2 + p[2]^2 + p[1] - t - h = (t, x, p) -> x[2]^2 + -2x[1]^2 + p[1]^2 - 2p[2]^2 + t - f₊g = (t, x, p) -> f(t, x, p) + g(t, x, p) - fg = (t, x, p) -> f(t, x, p) * g(t, x, p) - F = Hamiltonian(f; autonomous=false) - G = Hamiltonian(g; autonomous=false) - H = Hamiltonian(h; autonomous=false) - F₊G = Hamiltonian(f₊g; autonomous=false) - FG = Hamiltonian(fg; autonomous=false) - # - Test.@test Poisson(f, g; autonomous=false)(2, [1, 2], [2, 1]) == -28 - Test.@test Poisson(f, G)(2, [1, 2], [2, 1]) == -28 - Test.@test Poisson(F, g)(2, [1, 2], [2, 1]) == -28 - # - Test.@test Poisson(F, Hamiltonian((t, x, p) -> 42; autonomous=false))( - 2, [1, 2], [2, 1] - ) == 0 - Test.@test Poisson(F, G)(2, [1, 2], [2, 1]) == -28 - Test.@test Poisson(F, G)(2, [1, 2], [2, 1]) == -Poisson(G, F)(2, [1, 2], [2, 1]) # anticommutativity - Test.@test Poisson(F₊G, H)(2, [1, 2], [2, 1]) == - Poisson(F, H)(2, [1, 2], [2, 1]) + Poisson(G, H)(2, [1, 2], [2, 1]) # bilinearity 1 - Test.@test Poisson(H, F₊G)(2, [1, 2], [2, 1]) == - Poisson(H, F)(2, [1, 2], [2, 1]) + Poisson(H, G)(2, [1, 2], [2, 1]) # bilinearity 2 - Test.@test Poisson(FG, H)(2, [1, 2], [2, 1]) == - Poisson(F, H)(2, [1, 2], [2, 1]) * G(2, [1, 2], [2, 1]) + - F(2, [1, 2], [2, 1]) * Poisson(G, H)(2, [1, 2], [2, 1]) # Liebniz's rule - Test.@test Poisson(F, Poisson(G, H))(2, [1, 2], [2, 1]) + - Poisson(G, Poisson(H, F))(2, [1, 2], [2, 1]) + - Poisson(H, Poisson(F, G))(2, [1, 2], [2, 1]) == 0 # Jacobi identity - end - - @testset "autonomous nonfixed case" begin - f = (x, p, v) -> v[1] * x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] - g = (x, p, v) -> 3x[2]^2 + -x[1]^2 + p[2]^2 + p[1] - v[2] - h = (x, p, v) -> x[2]^2 + -2x[1]^2 + p[1]^2 - 2p[2]^2 + v[2] - f₊g = (x, p, v) -> f(x, p, v) + g(x, p, v) - fg = (x, p, v) -> f(x, p, v) * g(x, p, v) - F = Hamiltonian(f; variable=true) - G = Hamiltonian(g; variable=true) - H = Hamiltonian(h; variable=true) - F₊G = Hamiltonian(f₊g; variable=true) - FG = Hamiltonian(fg; variable=true) - Test.@test Poisson(F, Hamiltonian((x, p, v) -> 42; variable=true))( - [1, 2], [2, 1], [4, 4] - ) == 0 - Test.@test Poisson(F, G)([1, 2], [2, 1], [4, 4]) == -44 - Test.@test Poisson(F, G)([1, 2], [2, 1], [4, 4]) == - -Poisson(G, F)([1, 2], [2, 1], [4, 4]) # anticommutativity - Test.@test Poisson(F₊G, H)([1, 2], [2, 1], [4, 4]) == - Poisson(F, H)([1, 2], [2, 1], [4, 4]) + - Poisson(G, H)([1, 2], [2, 1], [4, 4]) # bilinearity 1 - Test.@test Poisson(H, F₊G)([1, 2], [2, 1], [4, 4]) == - Poisson(H, F)([1, 2], [2, 1], [4, 4]) + - Poisson(H, G)([1, 2], [2, 1], [4, 4]) # bilinearity 2 - Test.@test Poisson(FG, H)([1, 2], [2, 1], [4, 4]) == - Poisson(F, H)([1, 2], [2, 1], [4, 4]) * G([1, 2], [2, 1], [4, 4]) + - F([1, 2], [2, 1], [4, 4]) * Poisson(G, H)([1, 2], [2, 1], [4, 4]) # Liebniz's rule - Test.@test Poisson(F, Poisson(G, H))([1, 2], [2, 1], [4, 4]) + - Poisson(G, Poisson(H, F))([1, 2], [2, 1], [4, 4]) + - Poisson(H, Poisson(F, G))([1, 2], [2, 1], [4, 4]) == 0 # Jacobi identity - end - - @testset "nonautonomous nonfixed case" begin - f = (t, x, p, v) -> t * v[1] * x[2]^2 + 2x[1]^2 + p[1]^2 + v[2] - g = (t, x, p, v) -> 3x[2]^2 + -x[1]^2 + p[2]^2 + p[1] + t - v[2] - h = (t, x, p, v) -> x[2]^2 + -2x[1]^2 + p[1]^2 - 2p[2]^2 + t + v[2] - f₊g = (t, x, p, v) -> f(t, x, p, v) + g(t, x, p, v) - fg = (t, x, p, v) -> f(t, x, p, v) * g(t, x, p, v) - F = Hamiltonian(f; autonomous=false, variable=true) - G = Hamiltonian(g; autonomous=false, variable=true) - H = Hamiltonian(h; autonomous=false, variable=true) - F₊G = Hamiltonian(f₊g; autonomous=false, variable=true) - FG = Hamiltonian(fg; autonomous=false, variable=true) - Test.@test Poisson( - F, Hamiltonian((t, x, p, v) -> 42; autonomous=false, variable=true) - )( - 2, [1, 2], [2, 1], [4, 4] - ) == 0 - Test.@test Poisson(F, G)(2, [1, 2], [2, 1], [4, 4]) == -76 - Test.@test Poisson(F, G)(2, [1, 2], [2, 1], [4, 4]) == - -Poisson(G, F)(2, [1, 2], [2, 1], [4, 4]) # anticommutativity - Test.@test Poisson(F₊G, H)(2, [1, 2], [2, 1], [4, 4]) == - Poisson(F, H)(2, [1, 2], [2, 1], [4, 4]) + - Poisson(G, H)(2, [1, 2], [2, 1], [4, 4]) # bilinearity 1 - Test.@test Poisson(H, F₊G)(2, [1, 2], [2, 1], [4, 4]) == - Poisson(H, F)(2, [1, 2], [2, 1], [4, 4]) + - Poisson(H, G)(2, [1, 2], [2, 1], [4, 4]) # bilinearity 2 - Test.@test Poisson(FG, H)(2, [1, 2], [2, 1], [4, 4]) == - Poisson(F, H)(2, [1, 2], [2, 1], [4, 4]) * - G(2, [1, 2], [2, 1], [4, 4]) + - F(2, [1, 2], [2, 1], [4, 4]) * - Poisson(G, H)(2, [1, 2], [2, 1], [4, 4]) # Liebniz's rule - Test.@test Poisson(F, Poisson(G, H))(2, [1, 2], [2, 1], [4, 4]) + - Poisson(G, Poisson(H, F))(2, [1, 2], [2, 1], [4, 4]) + - Poisson(H, Poisson(F, G))(2, [1, 2], [2, 1], [4, 4]) == 0 # Jacobi identity - end + @testset "Variable dependence" begin + # H(x, p, v) + H(x, p, v) = v[1] * p[1] * x[1] + G(x, p, v) = x[1] + # {H, G} -> ∂H/∂p * ∂G/∂x ... + # ∂H/∂p = [v*x1, 0], ∂G/∂x = [1, 0] -> v*x1 + PB = CTFlows.Poisson(H, G; variable=true) + # v=2, x=[3, 4], p=[5, 6] -> 2*3 = 6 + @test PB([3.0, 4.0], [5.0, 6.0], [2.0]) ≈ 6.0 atol = 1e-6 + + # Macro with options + PB_macro = CTFlows.@Lie {H, G} variable = true + @test PB_macro([3.0, 4.0], [5.0, 6.0], [2.0]) ≈ 6.0 atol = 1e-6 end - @testset "poisson bracket of Lifts" begin - @testset "autonomous case" begin - f = x -> [x[1] + x[2]^2, x[1], 0] - g = x -> [0, x[2], x[1]^2 + 4 * x[2]] - F = Lift(f) - G = Lift(g) - F_ = (x, p) -> p' * f(x) - G_ = (x, p) -> p' * g(x) - Test.@test Poisson(F, G)([1, 2, 3], [4, 0, 4]) ≈ - Poisson(F_, G_)([1, 2, 3], [4, 0, 4]) atol = 1e-6 - Test.@test Poisson(F, G_)([1, 2, 3], [4, 0, 4]) ≈ - Poisson(F_, G)([1, 2, 3], [4, 0, 4]) atol = 1e-6 - end - - @testset "nonautonomous case" begin - f = (t, x) -> [t * x[1] + x[2]^2, x[1], 0] - g = (t, x) -> [0, x[2], t * x[1]^2 + 4 * x[2]] - F = Lift(f; autonomous=false, variable=false) - G = Lift(g; autonomous=false, variable=false) - F_ = (t, x, p) -> p' * f(t, x) - G_ = (t, x, p) -> p' * g(t, x) - Test.@test Poisson(F, G; autonomous=false, variable=false)( - 2, [1, 2, 3], [4, 0, 4] - ) ≈ Poisson(F_, G_; autonomous=false, variable=false)( - 2, [1, 2, 3], [4, 0, 4] - ) atol = 1e-6 - Test.@test Poisson(F, G_; autonomous=false, variable=false)( - 2, [1, 2, 3], [4, 0, 4] - ) ≈ Poisson(F_, G; autonomous=false, variable=false)( - 2, [1, 2, 3], [4, 0, 4] - ) atol = 1e-6 - end - - @testset "autonomous nonfixed case" begin - f = (x, v) -> [x[1] + v * x[2]^2, x[1], 0] - g = (x, v) -> [0, x[2], x[1]^2 + v * 4 * x[2]] - F = Lift(f; autonomous=true, variable=true) - G = Lift(g; autonomous=true, variable=true) - F_ = (x, p, v) -> p' * f(x, v) - G_ = (x, p, v) -> p' * g(x, v) - Test.@test Poisson(F, G; autonomous=true, variable=true)( - [1, 2, 3], [4, 0, 4], 1 - ) ≈ Poisson(F_, G_; autonomous=true, variable=true)( - [1, 2, 3], [4, 0, 4], 1 - ) atol = 1e-6 - Test.@test Poisson(F, G_; autonomous=true, variable=true)( - [1, 2, 3], [4, 0, 4], 1 - ) ≈ Poisson(F_, G; autonomous=true, variable=true)( - [1, 2, 3], [4, 0, 4], 1 - ) atol = 1e-6 - end - - @testset "nonautonomous nonfixed case" begin - f = (t, x, v) -> [t * x[1] + v * x[2]^2, x[1], 0] - g = (t, x, v) -> [0, x[2], t * x[1]^2 + v * 4 * x[2]] - F = Lift(f; autonomous=false, variable=true) - G = Lift(g; autonomous=false, variable=true) - F_ = (t, x, p, v) -> p' * f(t, x, v) - G_ = (t, x, p, v) -> p' * g(t, x, v) - Test.@test Poisson(F, G; autonomous=false, variable=true)( - 2, [1, 2, 3], [4, 0, 4], 1 - ) ≈ Poisson(F_, G_; autonomous=false, variable=true)( - 2, [1, 2, 3], [4, 0, 4], 1 - ) atol = 1e-6 - Test.@test Poisson(F, G_; autonomous=false, variable=true)( - 2, [1, 2, 3], [4, 0, 4], 1 - ) ≈ Poisson(F_, G; autonomous=false, variable=true)( - 2, [1, 2, 3], [4, 0, 4], 1 - ) atol = 1e-6 - end + @testset "Backend parameter" begin + using DifferentiationInterface + # Test that Poisson accepts backend parameter + H(x, p) = 0.5 * (p[1]^2 + p[2]^2) + G(x, p) = x[1] + backend = AutoForwardDiff() + PB = CTFlows.Poisson(H, G; backend=backend) + @test PB([1.0, 2.0], [3.0, 4.0]) ≈ 3.0 atol = 1e-6 + + # Test that ad accepts backend parameter + X(x) = [x[2], -x[1]] + f(x) = x[1]^2 + x[2]^2 + Lf = CTFlows.ad(X, f; backend=backend) + @test Lf([1.0, 2.0]) ≈ 0.0 atol = 1e-6 end - # macros + @testset "Nested Brackets (Jacobi Identity check)" begin + # [X, [Y, Z]] + [Y, [Z, X]] + [Z, [X, Y]] = 0 + X(x) = [0, x[3], -x[2]] + Y(x) = [-x[3], 0, x[1]] + Z(x) = [x[2], -x[1], 0] + + # We need to compose functions for nested brackets + # Let's verify expansion + XYZ = CTFlows.@Lie [X, [Y, Z]] + YZX = CTFlows.@Lie [Y, [Z, X]] + ZXY = CTFlows.@Lie [Z, [X, Y]] + + x0 = [1.0, 2.0, 3.0] + sum_jacobi = XYZ(x0) + YZX(x0) + ZXY(x0) + @test isapprox_vec(sum_jacobi, [0.0, 0.0, 0.0]) + end - @testset "lie macro" begin + @testset "Prefix System" begin + # Save current prefix + old_prefix = CTFlows.diffgeo_prefix() - # parameters - t = 1 - x = [1, 2, 3] - Γ = 2 - γ = 1 - δ = γ - Γ - v = 1 - - # autonomous - F0 = VectorField(x -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]) - F1 = VectorField(x -> [0, -x[3], x[2]]) - F2 = VectorField(x -> [x[3], 0, -x[1]]) - F01_ = Lie(F0, F1) - F011_ = Lie(F01_, F1) - F01__ = CTFlows.@Lie [F0, F1] - F011__ = CTFlows.@Lie [[F0, F1], F1] - - Test.@test F01_(x) ≈ F01__(x) atol = 1e-6 - Test.@test F011_(x) ≈ F011__(x) atol = 1e-6 - # - get_F0 = () -> F0 - F011___ = CTFlows.@Lie [[get_F0(), F1], F1] - Test.@test F011_(x) ≈ F011___(x) atol = 1e-6 - - # nonautonomous - F0 = VectorField( - (t, x) -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]; - autonomous=false, - variable=false, - ) - F1 = VectorField((t, x) -> [0, -x[3], x[2]]; autonomous=false, variable=false) - F2 = VectorField((t, x) -> [x[3], 0, -x[1]]; autonomous=false, variable=false) - F01_ = Lie(F0, F1) - F011_ = Lie(F01_, F1) - F01__ = CTFlows.@Lie [F0, F1] - F011__ = CTFlows.@Lie [[F0, F1], F1] - - Test.@test F01_(t, x) ≈ F01__(t, x) atol = 1e-6 - Test.@test F011_(t, x) ≈ F011__(t, x) atol = 1e-6 - # - get_F0 = () -> F0 - F011___ = CTFlows.@Lie [[get_F0(), F1], F1] - Test.@test F011_(t, x) ≈ F011___(t, x) atol = 1e-6 - - # autonomous nonfixed - F0 = VectorField( - (x, v) -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]; autonomous=true, variable=true - ) - F1 = VectorField((x, v) -> [0, -x[3], x[2]]; autonomous=true, variable=true) - F2 = VectorField((x, v) -> [x[3], 0, -x[1]]; autonomous=true, variable=true) - F01_ = Lie(F0, F1) - F011_ = Lie(F01_, F1) - F01__ = CTFlows.@Lie [F0, F1] - F011__ = CTFlows.@Lie [[F0, F1], F1] - - Test.@test F01_(x, v) ≈ F01__(x, v) atol = 1e-6 - Test.@test F011_(x, v) ≈ F011__(x, v) atol = 1e-6 - # - get_F0 = () -> F0 - F011___ = CTFlows.@Lie [[get_F0(), F1], F1] - Test.@test F011_(x, v) ≈ F011___(x, v) atol = 1e-6 - - # nonautonomous nonfixed - F0 = VectorField( - (t, x, v) -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]; - autonomous=false, - variable=true, - ) - F1 = VectorField((t, x, v) -> [0, -x[3], x[2]]; autonomous=false, variable=true) - F2 = VectorField((t, x, v) -> [x[3], 0, -x[1]]; autonomous=false, variable=true) - F01_ = Lie(F0, F1) - F011_ = Lie(F01_, F1) - F01__ = CTFlows.@Lie [F0, F1] - F011__ = CTFlows.@Lie [[F0, F1], F1] - - Test.@test F01_(t, x, v) ≈ F01__(t, x, v) atol = 1e-6 - Test.@test F011_(t, x, v) ≈ F011__(t, x, v) atol = 1e-6 - # - get_F0 = () -> F0 - F011___ = CTFlows.@Lie [[get_F0(), F1], F1] - Test.@test F011_(t, x, v) ≈ F011___(t, x, v) atol = 1e-6 + # Change prefix + CTFlows.diffgeo_prefix!(:MyModule) + @test CTFlows.diffgeo_prefix() == :MyModule + + # Check macro expansion uses new prefix (conceptually) + # We can't easily test macro expansion result availability unless MyModule exists + # But we tested the set/get mechanism + + # Restore prefix + CTFlows.diffgeo_prefix!(old_prefix) + @test CTFlows.diffgeo_prefix() == :CTFlows end - @testset "poisson macro" begin - # parameters - t = 1 - x = [1, 2, 3] - p = [1, 0, 7] - Γ = 2 - γ = 1 - δ = γ - Γ - v = 2 - - # autonomous - H0 = Hamiltonian((x, p) -> 0.5 * (x[1]^2 + x[2]^2 + p[1]^2)) - H1 = Hamiltonian((x, p) -> 0.5 * (x[1]^2 + x[2]^2 + p[2]^2)) - P01 = Poisson(H0, H1) - P011 = Poisson(P01, H1) - P01_ = CTFlows.@Lie {H0, H1} - P011_ = CTFlows.@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__ = CTFlows.@Lie {{get_H0(), H1}, H1} - Test.@test P011_(x, p) ≈ P011__(x, p) atol = 1e-6 - - # nonautonomous - H0 = Hamiltonian((t, x, p) -> 0.5 * (x[1]^2 + x[2]^2 + p[1]^2); autonomous=false) - H1 = Hamiltonian( - (t, x, p) -> 0.5 * (x[1]^2 + x[2]^2 + p[2]^2); autonomous=false, variable=false - ) - P01 = Poisson(H0, H1) - P011 = Poisson(P01, H1) - P01_ = CTFlows.@Lie {H0, H1} - P011_ = CTFlows.@Lie {{H0, H1}, H1} - Test.@test P01(t, x, p) ≈ P01_(t, x, p) atol = 1e-6 - Test.@test P011(t, x, p) ≈ P011_(t, x, p) atol = 1e-6 - get_H0 = () -> H0 - P011__ = CTFlows.@Lie {{get_H0(), H1}, H1} - Test.@test P011_(t, x, p) ≈ P011__(t, x, p) atol = 1e-6 - - # autonomous nonfixed - H0 = Hamiltonian((x, p, v) -> 0.5 * (x[1]^2 + x[2]^2 + p[1]^2 + v); variable=true) - H1 = Hamiltonian((x, p, v) -> 0.5 * (x[1]^2 + x[2]^2 + p[2]^2 + v); variable=true) - P01 = Poisson(H0, H1) - P011 = Poisson(P01, H1) - P01_ = CTFlows.@Lie {H0, H1} - P011_ = CTFlows.@Lie {{H0, H1}, H1} - 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__ = CTFlows.@Lie {{get_H0(), H1}, H1} - Test.@test P011_(x, p, v) ≈ P011__(x, p, v) atol = 1e-6 - - # nonautonomous nonfixed - H0 = Hamiltonian( - (t, x, p, v) -> 0.5 * (x[1]^2 + x[2]^2 + p[1]^2 + v); - autonomous=false, - variable=true, - ) - H1 = Hamiltonian( - (t, x, p, v) -> 0.5 * (x[1]^2 + x[2]^2 + p[2]^2 + v); - autonomous=false, - variable=true, - ) - P01 = Poisson(H0, H1) - P011 = Poisson(P01, H1) - P01_ = CTFlows.@Lie {H0, H1} - P011_ = CTFlows.@Lie {{H0, H1}, H1} - 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__ = CTFlows.@Lie {{get_H0(), H1}, H1} - Test.@test P011_(t, x, p, v) ≈ P011__(t, x, p, v) atol = 1e-6 + @testset "∂ₜ time derivative" begin + f(t, x) = t^2 * x[1] + df = CTFlows.∂ₜ(f) + # ∂f/∂t = 2t * x1 + @test df(3.0, [2.0]) ≈ 2 * 3.0 * 2.0 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_ = CTFlows.@Lie {H0, H1} - P011_ = CTFlows.@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__ = CTFlows.@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 = CTFlows.@Lie {H0, H1}(t, x, p) autonomous = false - P01_ = CTFlows.@Lie {H0, H1} autonomous = false - P011_ = CTFlows.@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__ = CTFlows.@Lie {{get_H0(), H1}, H1} autonomous = false variable = false - 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_ = CTFlows.@Lie {H0, H1} variable = true - P011_ = CTFlows.@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__ = CTFlows.@Lie {{get_H0(), H1}, H1} autonomous = true variable = true - 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_ = CTFlows.@Lie {H0, H1} autonomous = false variable = true - P011_ = CTFlows.@Lie {{H0, H1}, H1} autonomous = false variable = true - 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__ = CTFlows.@Lie {{get_H0(), H1}, H1} autonomous = false variable = true - Test.@test P011_(t, x, p, v) ≈ P011__(t, x, p, v) atol = 1e-6 + # ============================================================================ + # PRIORITY 1: Type-Based API Tests + # ============================================================================ + + @testset "Type-Based API - ad()" begin + # Test that explicit type dispatch works correctly + X(x) = [x[2], -x[1]] + f(x) = x[1]^2 + x[2]^2 + + # Autonomous, Fixed + Lf_typed = CTFlows.ad(X, f, CTFlows.Autonomous, CTFlows.Fixed) + Lf_kwargs = CTFlows.ad(X, f; autonomous=true, variable=false) + @test Lf_typed([1.0, 2.0]) ≈ Lf_kwargs([1.0, 2.0]) atol = 1e-6 + + # NonAutonomous, Fixed + X_na(t, x) = [t * x[2], -x[1]] + f_na(t, x) = t + x[1]^2 + Lf_na_typed = CTFlows.ad(X_na, f_na, CTFlows.NonAutonomous, CTFlows.Fixed) + Lf_na_kwargs = CTFlows.ad(X_na, f_na; autonomous=false, variable=false) + @test Lf_na_typed(2.0, [1.0, 2.0]) ≈ Lf_na_kwargs(2.0, [1.0, 2.0]) atol = 1e-6 + + # Autonomous, NonFixed (with variable) + X_v(x, v) = [v[1] * x[2], -x[1]] + f_v(x, v) = v[1] * x[1]^2 + x[2]^2 + Lf_v_typed = CTFlows.ad(X_v, f_v, CTFlows.Autonomous, CTFlows.NonFixed) + Lf_v_kwargs = CTFlows.ad(X_v, f_v; autonomous=true, variable=true) + @test Lf_v_typed([1.0, 2.0], [2.0]) ≈ Lf_v_kwargs([1.0, 2.0], [2.0]) atol = 1e-6 + + # NonAutonomous, NonFixed + X_nav(t, x, v) = [t * v[1] * x[2], -x[1]] + f_nav(t, x, v) = t + v[1] * x[1]^2 + Lf_nav_typed = CTFlows.ad(X_nav, f_nav, CTFlows.NonAutonomous, CTFlows.NonFixed) + Lf_nav_kwargs = CTFlows.ad(X_nav, f_nav; autonomous=false, variable=true) + @test Lf_nav_typed(2.0, [1.0, 2.0], [2.0]) ≈ Lf_nav_kwargs(2.0, [1.0, 2.0], [2.0]) atol = 1e-6 + end + + @testset "Type-Based API - Lift()" begin + # Autonomous, Fixed + f(x) = [x[1]^2, x[2]] + H_typed = CTFlows.Lift(f, CTFlows.Autonomous, CTFlows.Fixed) + H_kwargs = CTFlows.Lift(f; autonomous=true, variable=false) + @test H_typed([1.0, 2.0], [3.0, 4.0]) ≈ H_kwargs([1.0, 2.0], [3.0, 4.0]) atol = 1e-6 + + # NonAutonomous, Fixed + g(t, x) = [t * x[1], x[2]] + Hg_typed = CTFlows.Lift(g, CTFlows.NonAutonomous, CTFlows.Fixed) + Hg_kwargs = CTFlows.Lift(g; autonomous=false, variable=false) + @test Hg_typed(2.0, [1.0, 2.0], [3.0, 4.0]) ≈ Hg_kwargs(2.0, [1.0, 2.0], [3.0, 4.0]) atol = 1e-6 + + # Autonomous, NonFixed + h(x, v) = [v[1] * x[1], x[2]] + Hh_typed = CTFlows.Lift(h, CTFlows.Autonomous, CTFlows.NonFixed) + Hh_kwargs = CTFlows.Lift(h; autonomous=true, variable=true) + @test Hh_typed([1.0, 2.0], [3.0, 4.0], [2.0]) ≈ Hh_kwargs([1.0, 2.0], [3.0, 4.0], [2.0]) atol = 1e-6 + end + + @testset "Type-Based API - Poisson()" begin + # Autonomous, Fixed + H(x, p) = 0.5 * (p[1]^2 + p[2]^2) + G(x, p) = x[1] + PB_typed = CTFlows.Poisson(H, G, CTFlows.Autonomous, CTFlows.Fixed) + PB_kwargs = CTFlows.Poisson(H, G; autonomous=true, variable=false) + @test PB_typed([1.0, 2.0], [3.0, 4.0]) ≈ PB_kwargs([1.0, 2.0], [3.0, 4.0]) atol = 1e-6 + + # NonAutonomous, Fixed + H_na(t, x, p) = t + 0.5 * (p[1]^2 + p[2]^2) + G_na(t, x, p) = x[1] + PB_na_typed = CTFlows.Poisson(H_na, G_na, CTFlows.NonAutonomous, CTFlows.Fixed) + PB_na_kwargs = CTFlows.Poisson(H_na, G_na; autonomous=false, variable=false) + @test PB_na_typed(1.0, [1.0, 2.0], [3.0, 4.0]) ≈ PB_na_kwargs(1.0, [1.0, 2.0], [3.0, 4.0]) atol = 1e-6 + end + + # ============================================================================ + # PRIORITY 2: Edge Cases Tests + # ============================================================================ + + @testset "Edge Cases - 1D (Scalar)" begin + # ad() with scalars + X_1d(x) = -x + f_1d(x) = x^2 + Lf_1d = CTFlows.ad(X_1d, f_1d) + # L_X(f) = df/dx * X = 2x * (-x) = -2x^2 + # At x=3: -2*9 = -18 + @test Lf_1d(3.0) ≈ -18.0 atol = 1e-6 + + # Poisson() with scalars (already tested above, but explicit) + H_1d(x, p) = 0.5 * (p^2 + x^2) + G_1d(x, p) = x + PB_1d = CTFlows.Poisson(H_1d, G_1d) + @test PB_1d(1.0, 3.0) ≈ 3.0 atol = 1e-6 + + # Lift() with scalar + f_1d_vec(x) = -2 * x # returns scalar (treated as 1D vector) + H_1d = CTFlows.Lift(f_1d_vec) + @test H_1d(2.0, 3.0) ≈ -12.0 atol = 1e-6 + end + + @testset "Edge Cases - High Dimensional (3D, 4D)" begin + # 3D Lie bracket + X_3d(x) = [0, x[3], -x[2]] + Y_3d(x) = [-x[3], 0, x[1]] + XY_3d = CTFlows.ad(X_3d, Y_3d) + # Known result for rotational vector fields + @test length(XY_3d([1.0, 2.0, 3.0])) == 3 + + # 4D Poisson bracket + H_4d(x, p) = sum(p .^ 2) + sum(x .^ 2) + G_4d(x, p) = x[1] * p[1] + PB_4d = CTFlows.Poisson(H_4d, G_4d) + result_4d = PB_4d([1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0]) + @test result_4d isa Number # Just check it computes + end + + @testset "Edge Cases - Trivial Cases" begin + # Zero vector field + X_zero(x) = [0.0, 0.0] + f(x) = x[1]^2 + x[2]^2 + Lf_zero = CTFlows.ad(X_zero, f) + @test Lf_zero([1.0, 2.0]) ≈ 0.0 atol = 1e-10 + + # Constant function (gradient = 0) + X(x) = [x[2], -x[1]] + f_const(x) = 42.0 + Lf_const = CTFlows.ad(X, f_const) + @test Lf_const([1.0, 2.0]) ≈ 0.0 atol = 1e-10 + + # Zero Hamiltonian + H_zero(x, p) = 0.0 + G(x, p) = x[1] + PB_zero = CTFlows.Poisson(H_zero, G) + @test PB_zero([1.0, 2.0], [3.0, 4.0]) ≈ 0.0 atol = 1e-10 + end + + @testset "Edge Cases - Commuting Vector Fields" begin + # [X, Y] = 0 for commuting fields + X(x) = [1.0, 0.0] # Constant field in x direction + Y(x) = [0.0, 1.0] # Constant field in y direction + XY = CTFlows.ad(X, Y) + # Commuting vector fields should have zero Lie bracket + @test isapprox_vec(XY([1.0, 2.0]), [0.0, 0.0]) + end + + # ============================================================================ + # PACKAGE 1: Mathematical Properties (from backup analysis) + # ============================================================================ + + @testset "Poisson - Mathematical Properties" begin + # Setup test functions + f(x, p) = x[2]^2 + 2x[1]^2 + p[1]^2 + g(x, p) = 3x[2]^2 + -x[1]^2 + p[2]^2 + p[1] + h(x, p) = x[2]^2 + -2x[1]^2 + p[1]^2 - 2p[2]^2 + f_plus_g(x, p) = f(x, p) + g(x, p) + f_times_g(x, p) = f(x, p) * g(x, p) + const_42(x, p) = 42.0 + + x_test = [1, 2] + p_test = [2, 1] + + # Property 1: Constant function has zero Poisson bracket + @test CTFlows.Poisson(f, const_42)(x_test, p_test) ≈ 0.0 atol = 1e-10 + + # Property 2: Anticommutativity - {F, G} = -{G, F} + PB_fg = CTFlows.Poisson(f, g)(x_test, p_test) + PB_gf = CTFlows.Poisson(g, f)(x_test, p_test) + @test PB_fg ≈ -PB_gf atol = 1e-10 + + # Property 3: Bilinearity (left) - {F+G, H} = {F, H} + {G, H} + PB_fpg_h = CTFlows.Poisson(f_plus_g, h)(x_test, p_test) + PB_f_h = CTFlows.Poisson(f, h)(x_test, p_test) + PB_g_h = CTFlows.Poisson(g, h)(x_test, p_test) + @test PB_fpg_h ≈ PB_f_h + PB_g_h atol = 1e-10 + + # Property 4: Bilinearity (right) - {H, F+G} = {H, F} + {H, G} + PB_h_fpg = CTFlows.Poisson(h, f_plus_g)(x_test, p_test) + PB_h_f = CTFlows.Poisson(h, f)(x_test, p_test) + PB_h_g = CTFlows.Poisson(h, g)(x_test, p_test) + @test PB_h_fpg ≈ PB_h_f + PB_h_g atol = 1e-10 + + # Property 5: Leibniz's rule - {FG, H} = {F, H}·G + F·{G, H} + PB_ftg_h = CTFlows.Poisson(f_times_g, h)(x_test, p_test) + leibniz_rhs = PB_f_h * g(x_test, p_test) + f(x_test, p_test) * PB_g_h + @test PB_ftg_h ≈ leibniz_rhs atol = 1e-10 + + # Property 6: Jacobi identity - {F, {G, H}} + {G, {H, F}} + {H, {F, G}} = 0 + PB_gh = CTFlows.Poisson(g, h) + PB_hf = CTFlows.Poisson(h, f) + PB_fg_func = CTFlows.Poisson(f, g) + + PB_f_gh = CTFlows.Poisson(f, PB_gh)(x_test, p_test) + PB_g_hf = CTFlows.Poisson(g, PB_hf)(x_test, p_test) + PB_h_fg = CTFlows.Poisson(h, PB_fg_func)(x_test, p_test) + + jacobi_sum = PB_f_gh + PB_g_hf + PB_h_fg + @test abs(jacobi_sum) < 1e-10 end - @testset "lie and poisson macros operation (sum,diff,...)" begin - # parameters - t = 1 - x = [1, 2, 3] - p = [1, 0, 7] - Γ = 2 - γ = 1 + @testset "Poisson of Lifts - Composition" begin + # Test that Poisson(Lift(f), Lift(g)) gives correct results + f(x) = [x[1] + x[2]^2, x[1], 0] + g(x) = [0, x[2], x[1]^2 + 4 * x[2]] + + # Create Lifts + F = CTFlows.Lift(f) + G = CTFlows.Lift(g) + + # Explicit Hamiltonians + F_explicit(x, p) = p' * f(x) + G_explicit(x, p) = p' * g(x) + + x_test = [1.0, 2.0, 3.0] + p_test = [4.0, 0.0, 4.0] + + # Poisson of Lifts should equal Poisson of explicit Hamiltonians + @test CTFlows.Poisson(F, G)(x_test, p_test) ≈ + CTFlows.Poisson(F_explicit, G_explicit)(x_test, p_test) atol = 1e-6 + + # Mixed case: Lift + explicit + @test CTFlows.Poisson(F, G_explicit)(x_test, p_test) ≈ + CTFlows.Poisson(F_explicit, G)(x_test, p_test) atol = 1e-6 + + # Non-autonomous case + f_na(t, x) = [t * x[1] + x[2]^2, x[1], 0] + g_na(t, x) = [0, x[2], t * x[1]^2 + 4 * x[2]] + + F_na = CTFlows.Lift(f_na; autonomous=false) + G_na = CTFlows.Lift(g_na; autonomous=false) + F_na_explicit(t, x, p) = p' * f_na(t, x) + G_na_explicit(t, x, p) = p' * g_na(t, x) + + t_test = 2.0 + @test CTFlows.Poisson(F_na, G_na; autonomous=false)(t_test, x_test, p_test) ≈ + CTFlows.Poisson(F_na_explicit, G_na_explicit; autonomous=false)(t_test, x_test, p_test) atol = 1e-6 + end + + # ============================================================================ + # PACKAGE 2: Physical Examples & Theoretical Validation + # ============================================================================ + + @testset "MRI Example - Bloch Equations" begin + # Physical constants for magnetic resonance imaging + Γ = 2.0 # Relaxation rate + γ = 1.0 # Gyromagnetic ratio δ = γ - Γ - v = 1 - - # lie - # autonomous - F0 = VectorField(x -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]) - F1 = VectorField(x -> [0, -x[3], x[2]]) - F2 = VectorField(x -> [x[3], 0, -x[1]]) - Test.@test CTFlows.@Lie [F0, F1](x) + 4 * [F1, F2](x) == [8, -8, -2] - Test.@test CTFlows.@Lie [F0, F1](x) - [F1, F2](x) == [-2, -3, -2] - Test.@test CTFlows.@Lie [F0, F1](x) .* [F1, F2](x) == [0, 4, 0] - Test.@test CTFlows.@Lie [1, 1, 1] + ([[F0, F1], F1](x) + [F1, F2](x) + [1, 1, 1]) == - [4, 5, -5] - - # nonautonomous nonfixed - F0 = VectorField( - (t, x, v) -> [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])]; - autonomous=false, - variable=true, - ) - F1 = VectorField((t, x, v) -> [0, -x[3], x[2]]; autonomous=false, variable=true) - F2 = VectorField((t, x, v) -> [x[3], 0, -x[1]]; autonomous=false, variable=true) - Test.@test CTFlows.@Lie [F0, F1](t, x, v) + 4 * [F1, F2](t, x, v) == [8, -8, -2] - Test.@test CTFlows.@Lie [F0, F1](t, x, v) - [F1, F2](t, x, v) == [-2, -3, -2] - Test.@test CTFlows.@Lie [F0, F1](t, x, v) .* [F1, F2](t, x, v) == [0, 4, 0] - Test.@test CTFlows.@Lie [1, 1, 1] + - ([[F0, F1], F1](t, x, v) + [F1, F2](t, x, v) + [1, 1, 1]) == - [4, 5, -5] - - # poisson - # autonomous - H0 = Hamiltonian((x, p) -> 0.5 * (2x[1]^2 + x[2]^2 + p[1]^2)) - H1 = Hamiltonian((x, p) -> 0.5 * (3x[1]^2 + x[2]^2 + p[2]^2)) - H2 = Hamiltonian((x, p) -> 0.5 * (4x[1]^2 + x[2]^2 + p[1]^3 + p[2]^2)) - Test.@test CTFlows.@Lie {H0, H1}(x, p) + 4 * {H1, H2}(x, p) == -15 - Test.@test CTFlows.@Lie {H0, H1}(x, p) - {H1, H2}(x, p) == 7.5 - Test.@test CTFlows.@Lie {H0, H1}(x, p) * {H1, H2}(x, p) == -13.5 - Test.@test CTFlows.@Lie 4 + ({{H0, H1}, H1}(x, p) + -2 * {H1, H2}(x, p) + 21) == 39 - - # nonautonomous nonfixed - H0 = Hamiltonian( - (t, x, p, v) -> 0.5 * (2x[1]^2 + x[2]^2 + p[1]^2); - autonomous=false, - variable=true, - ) - H1 = Hamiltonian( - (t, x, p, v) -> 0.5 * (3x[1]^2 + x[2]^2 + p[2]^2); - autonomous=false, - variable=true, - ) - H2 = Hamiltonian( - (t, x, p, v) -> 0.5 * (4x[1]^2 + x[2]^2 + p[1]^3 + p[2]^2); - autonomous=false, - variable=true, - ) - Test.@test CTFlows.@Lie {H0, H1}(t, x, p, v) + 4 * {H1, H2}(t, x, p, v) == -15 - Test.@test CTFlows.@Lie {H0, H1}(t, x, p, v) - {H1, H2}(t, x, p, v) == 7.5 - Test.@test CTFlows.@Lie {H0, H1}(t, x, p, v) * {H1, H2}(t, x, p, v) == -13.5 - Test.@test CTFlows.@Lie 4 + ( - {{H0, H1}, H1}(t, x, p, v) + -2 * {H1, H2}(t, x, p, v) + 21 - ) == 39 + + # Bloch equation vector fields + F0(x) = [-Γ * x[1], -Γ * x[2], γ * (1 - x[3])] + F1(x) = [0.0, -x[3], x[2]] + F2(x) = [x[3], 0.0, -x[1]] + + # Compute Lie brackets + F01 = CTFlows.ad(F0, F1) + F02 = CTFlows.ad(F0, F2) + F12 = CTFlows.ad(F1, F2) + + x = [1.0, 2.0, 3.0] + + # Verify known analytical results for Bloch equations + # WHY: These are well-known results in MRI physics + @test F01(x) ≈ -[0.0, γ - δ * x[3], -δ * x[2]] atol = 1e-6 + @test F02(x) ≈ -[-γ + δ * x[3], 0.0, δ * x[1]] atol = 1e-6 + @test F12(x) ≈ -[-x[2], x[1], 0.0] atol = 1e-6 end -end # test_differential_geometry + + @testset "Lie Bracket - Intrinsic Definition" begin + # Verify intrinsic definition: [X, Y]·f = X·(Y·f) - Y·(X·f) + # WHY: This is the fundamental commutator property of Lie brackets + + X(x) = [x[2]^2, -2x[1] * x[2]] + Y(x) = [x[1] * (1 + x[2]), 3x[2]^3] + f(x) = x[1]^4 + 2x[2]^3 + + x_test = [1.0, 2.0] + + # Method 1: Direct computation of [X,Y]·f + XY = CTFlows.ad(X, Y) + XY_dot_f = CTFlows.ad(XY, f) + result_direct = XY_dot_f(x_test) + + # Method 2: Commutator of directional derivatives X·(Y·f) - Y·(X·f) + Y_dot_f = CTFlows.ad(Y, f) + X_dot_f = CTFlows.ad(X, f) + X_dot_Yf = CTFlows.ad(X, Y_dot_f) + Y_dot_Xf = CTFlows.ad(Y, X_dot_f) + result_commutator = X_dot_Yf(x_test) - Y_dot_Xf(x_test) + + # Both methods should give the same result + @test result_direct ≈ result_commutator atol = 1e-6 + end + +end