[Dev] Refactoring

[ocots]

CTFlows.jl - Architecture and Use Cases (V3)

This document describes the architecture and use cases for CTFlows.jl version 3.

---

1. Overview

CTFlows.jl provides a unified interface for computing flows (time evolution) of dynamical systems arising from:

• Optimal Control Problems (OCP): Hamiltonian flows with various control specifications
• General ODEs: Systems defined via SciMLBase.ODEFunction

The package emphasizes:

• Simplicity: Minimal API surface, no wrappers or descriptions
• Extensibility: Backend-agnostic core with SciML integration via extensions
• Type Safety: Compile-time guarantees for flow concatenation
• Performance: Lightweight dependencies, efficient multi-phase flows

---

2. Use Cases

2.1 Flow from an Optimal Control Problem (OCP)

The primary use case is computing flows from an OCP with different control specifications.

UC1.1: Dynamic Feedback (Default)

ocp = @def ... # Define OCP
u(x, p) = p[2]  # Control depends on state and costate
f = Flow(ocp, u, :dynamic_feedback)
# or simply: f = Flow(ocp, u)  # defaults to :dynamic_feedback

# Point evaluation
xf, pf = f(t0, x0, p0, tf, v)

# Complete solution
sol = f((t0, tf), x0, p0, v)

UC1.2: State Feedback

u(x) = -0.5 * x  # Control depends only on state
f = Flow(ocp, u, :state_feedback)

# Usage
xf, pf = f(t0, x0, p0, tf, v)

UC1.3: Open Loop

u(t) = sin(t)  # Control depends only on time
f = Flow(ocp, u, :open_loop)

# Usage
xf, pf = f(t0, x0, p0, tf, v)

UC1.4: Regular Case with Constraints

# With state constraints and multipliers
f = Flow(ocp, u, g, μ, :dynamic_feedback)

# Usage
xf, pf = f(t0, x0, p0, tf, v)

UC1.5: Regular Case via Implicit Function Theorem

# OCP only, returns flow and ∂H/∂u
f, ∂H∂u = Flow(ocp)

# Usage: requires initial control guess
u0 = [0.1]
xf, pf, uf = f(t0, x0, p0, u0, tf, v)

---

2.2 Flow from an ODEFunction

For general ODEs already defined using SciML's ODEFunction. We support both out-of-place and in-place signatures.

using SciMLBase

# Out-of-place
of1 = ODEFunction((u, p, t) -> -u)
f1 = Flow(of1)

# In-place
of2 = ODEFunction((du, u, p, t) -> du .= -u)
f2 = Flow(of2)

# Usage
xf1 = f1(t0, x0, tf)
xf2 = f2(t0, x0, tf)

---

2.3 Options Management

Options are lightweight and can be defined at creation or overridden at call time.

# Creation with options
f = Flow(ocp, u, :state_feedback, alg=Tsit5(), abstol=1e-8)

# Call with overrides
sol = f((t0, tf), x0, p0, v; abstol=1e-10)

---

2.4 Flow Concatenation

Type Constraint

Flows can only be concatenated if they are of the same strict type:

• Two OCP flows must have the same control mode (both :state_feedback, both :dynamic_feedback, etc.)
• Two ODEFunction flows must come from the same ODEFunction

This constraint is enforced via type parameters that track the flow's origin.

Concatenation Syntax

# Valid: same control mode
f1 = Flow(ocp, u1, :state_feedback)
f2 = Flow(ocp, u2, :state_feedback)
f = f1 * (t1, f2)  # ✓ OK

# Invalid: different control modes
f1 = Flow(ocp, u1, :state_feedback)
f2 = Flow(ocp, u2, :dynamic_feedback)
f = f1 * (t1, f2)  # ✗ Type error

# With jumps
f = f1 * (t1, Δx, Δp, f2)  # Concatenation with state/costate jumps

Outputs

• Point evaluation: Returns the final state (and costate if applicable)
• Complete solution: Returns a merged solution object

xf, pf = f(t0, x0, p0, tf, v)      # Point
sol    = f((t0, tf), x0, p0, v)     # Merged solution

---

2.5 Autonomy and Variable Handling

Automatic Detection from OCP

The flow call signature is automatically determined from the OCP properties:

• Autonomy: Detected via CTModels.is_autonomous(ocp)
• Variable: Detected via CTModels.variable_dimension(ocp) > 0

The control function u must respect the OCP's autonomy:

• Autonomous OCP: u(x), u(x, p), or u(x, p, v)
• Non-Autonomous OCP: u(t, x), u(t, x, p), or u(t, x, p, v)

Variable Argument v

The variable v in flow calls depends on the OCP configuration:

Case 1: No variable (variable_dimension(ocp) == 0)

ocp = @def ... # No free times, no external variables
f = Flow(ocp, u, :state_feedback)

# Call without v
xf, pf = f(t0, x0, p0, tf)  # v is omitted

Case 2: Variable is tf (has_free_final_time(ocp) && variable_dimension(ocp) == 1)

ocp = @def ... # Free final time
f = Flow(ocp, u, :state_feedback)

# tf is already provided, no need for v
xf, pf = f(t0, x0, p0, tf)  # v = tf is inferred automatically

Case 3: Variable is t0 (has_free_initial_time(ocp) && variable_dimension(ocp) == 1)

ocp = @def ... # Free initial time
f = Flow(ocp, u, :state_feedback)

# t0 is already provided, no need for v
xf, pf = f(t0, x0, p0, tf)  # v = t0 is inferred automatically

Case 4: Variable is [t0, tf] (has_free_initial_time(ocp) && has_free_final_time(ocp) && variable_dimension(ocp) == 2)

ocp = @def ... # Both times free
f = Flow(ocp, u, :state_feedback)

# Both times are provided, no need for v
xf, pf = f(t0, x0, p0, tf)  # v = [t0, tf] is inferred automatically

Case 5: External variable (other cases)

ocp = @def ... # Has external parameter v
f = Flow(ocp, u, :state_feedback)

# Must provide v explicitly
v = [1.0, 2.0]
xf, pf = f(t0, x0, p0, tf, v)

Overriding Autonomy

You can override the OCP's autonomy for the control function using named arguments:

# Autonomous OCP, but time-dependent control
ocp = @def ... # Autonomous
u(t, x) = sin(t) * x  # Time-dependent control

f = Flow(ocp, u, :state_feedback, autonomous=false)  # Override

# Now the control is treated as non-autonomous
xf, pf = f(t0, x0, p0, tf)

Similarly, you can override the variable dependency:

# OCP without variables, but control depends on external parameter
ocp = @def ... # No variables
u(x, v) = v[1] * x  # Depends on external v

f = Flow(ocp, u, :state_feedback, variable=true)  # Override

# Now must provide v
v = [2.0]
xf, pf = f(t0, x0, p0, tf, v)

---

2.6 Differential Geometry Tools

CTFlows.jl provides a comprehensive library of differential geometry primitives for optimal control. These tools enable direct implementation of Pontryagin's Maximum Principle formulas without manual computation of partial derivatives.

Design Philosophy: No Type Wrappers

All differential geometry tools work with pure Function objects — there are no type wrappers like Hamiltonian, VectorField, or HamiltonianLift. This design choice provides:

1. Simplicity: No type hierarchy to learn or manage
2. Flexibility: All functions are composable Function objects
3. Clarity: Explicit function calls instead of operator overloading on custom types
4. Performance: No runtime type dispatch overhead

The only "types" are mathematical conventions:

• A vector field is a Function that maps x → ℝⁿ
• A Hamiltonian is a Function that maps (x, p) → ℝ
• A scalar function is a Function that maps x → ℝ

Available Tools

Tool	Notation	Signature	Description	Goddard
Hamiltonian Lift	Lift(f)	Function → Function	Lift f: (x, p) → p' * f(x)	✅
Hamiltonian Lift	Lift(f; autonomous, variable)	Function → Function	Lift with autonomy/variable options	✅
Lie Derivative	Lie(X, f; autonomous, variable)	(Function, Function) → Function	(X⋅f)(x) = ∇f(x)' * X(x)	✅
Time Derivative	∂ₜ(f)	Function → Function	(t, args...) → ∂f/∂t	❌
Lie Bracket	ad(X, Y; autonomous, variable)	(Function, Function) → Function	[X,Y](x) = J_Y(x)·X(x) - J_X(x)·Y(x)	❌
Poisson Bracket	Poisson(H, G; autonomous, variable)	(Function, Function) → Function	{H,G}(x,p) = ∇ₚH'·∇ₓG - ∇ₓH'·∇ₚG	❌
Lie Bracket Macro	@Lie [X, Y]	Macro	Lie bracket (expands to ad(X, Y))	❌
Poisson Bracket Macro	@Lie {H, G}	Macro	Poisson bracket (mathematical notation)	✅
Macro with Options	@Lie {H, G} autonomous=false	Macro	Macro with autonomy/variable options	❌

Example: Goddard Problem

The Goddard tutorial demonstrates practical usage of these tools for computing singular and boundary arc controls.

1. Hamiltonian Lift

H0 = Lift(F0)    # H0(x, p) = p' * F0(x)
H1 = Lift(F1)    # H1(x, p) = p' * F1(x)

2. Poisson Brackets for Singular Control

H01  = @Lie {H0, H1}        # {H0, H1}
H001 = @Lie {H0, H01}       # {H0, {H0, H1}}
H101 = @Lie {H1, H01}       # {H1, {H0, H1}}

# Singular control (minimal order)
us(x, p) = -H001(x, p) / H101(x, p)

3. Lie Derivative for Boundary Control

g(x) = vmax - x[2]  # State constraint

# Boundary control maintaining g(x) = 0
ub(x) = -Lie(F0, g, autonomous=true)(x) / Lie(F1, g, autonomous=true)(x)

4. Multiplier for State Constraint

μ(x, p) = H01(x, p) / Lie(F1, g, autonomous=true)(x)

Mathematical Formulas

Hamiltonian Lift

$$H_f(x, p) := p^\top f(x)$$

Lie Derivative

$$(X \cdot f)(x) := \nabla f(x)^\top X(x)$$

Lie Bracket

$$[X, Y](x) := J_Y(x) \cdot X(x) - J_X(x) \cdot Y(x)$$

Poisson Bracket

$$\{H, G\}(x, p) := \nabla_p H^\top \nabla_x G - \nabla_x H^\top \nabla_p G$$

where $J_Y$ is the Jacobian matrix of $Y$.

Implementation Details

No Type Dispatch: All operators work directly on Function objects:

# Lie derivative (explicit autonomous/variable arguments)
function Lie(X::Function, f::Function; autonomous=true, variable=false)::Function
    if autonomous
        return (x, args...) -> ctgradient(y -> f(y, args...), x)' * X(x, args...)
    else
        return (t, x, args...) -> ctgradient(y -> f(t, y, args...), x)' * X(t, x, args...)
    end
end

# Poisson bracket
function Poisson(H::Function, G::Function; autonomous=true, variable=false)::Function
    # Returns a function computing {H, G}
    # Implementation uses automatic differentiation (ForwardDiff.jl)
end

Automatic Differentiation: All derivatives are computed using ForwardDiff.jl via the ctgradient and ctjacobian utilities.

Macro Magic: The @Lie macro provides mathematical notation while dispatching to the appropriate function:

@Lie [X, Y]    # Expands to: ad(X, Y)
@Lie {H, G}    # Expands to: Poisson(H, G)

---

3. Internal Machinery

This section describes the internal architecture using Option 1: Implicit Dispatch on Algorithm Type.

3.1 Core Design Principles

1. Minimal Dependency: SciMLBase.jl extension only (triggered by using SciMLBase)
2. Stub Functions: _ode_problem and _solve are stubs in src/ that error with helpful messages
3. Extension Activation: Loading SciMLBase (via OrdinaryDiffEq, SimpleDiffEq, or DifferentialEquations) activates the extension
4. Multi-Phase Flows: Concatenated flows are represented as multi-phase systems with automatic phase transitions

3.2 Extension Loading Strategy

Stubs in src/CTFlows.jl

# src/stubs.jl
function _ode_problem(system, tspan, z0, params)
    @info """
    CTFlows.jl requires SciMLBase to create ODE problems.
    Please load a SciML solver package:
      using OrdinaryDiffEq  # Full solver suite
      using SimpleDiffEq    # Lightweight solvers
      using DifferentialEquations  # Complete ecosystem
    """
    error("SciMLBase extension not loaded. See message above.")
end

function _solve(problem, alg, options)
    @info """
    CTFlows.jl requires SciMLBase to solve ODE problems.
    Please load a SciML solver package:
      using OrdinaryDiffEq
      using SimpleDiffEq
      using DifferentialEquations
    """
    error("SciMLBase extension not loaded. See message above.")
end

Extension Override in ext/CTFlowsSciMLBaseExt.jl

# ext/CTFlowsSciMLBaseExt.jl
module CTFlowsSciMLBaseExt

using CTFlows
using SciMLBase

# Override stubs
function CTFlows._ode_problem(system::CTFlows.AbstractSystem, tspan, z0, params)
    rhs! = CTFlows.build_rhs(system)
    return ODEProblem(rhs!, z0, tspan, params)
end

function CTFlows._solve(problem::ODEProblem, alg, options)
    return solve(problem, alg; options...)
end

end

---

3.3 Multi-Phase Flow Concatenation

When flows are concatenated, a MultiPhaseSystem is created that manages phase transitions.

Internal Representation

# src/concatenation.jl
struct PhaseTransition
    time::Float64           # Transition time
    jump::Tuple{Vector, Vector}  # (Δx, Δp) jumps
end

struct MultiPhaseSystem{S<:AbstractSystem} <: AbstractSystem
    phases::Vector{S}       # List of subsystems
    transitions::Vector{PhaseTransition}  # Phase transition points
end

Point Evaluation Logic

For point evaluation xf, pf = f(t0, x0, p0, tf, v):

function evaluate_point(mps::MultiPhaseSystem, t0, z0, tf, v)
    # 1. Determine starting phase
    phase_idx = find_phase(mps, t0)
    
    # 2. Integrate through phases
    z_current = z0
    t_current = t0
    
    for k in phase_idx:length(mps.phases)
        # Determine integration endpoint for this phase
        t_end = k < length(mps.phases) ? mps.transitions[k].time : tf
        t_end = min(t_end, tf)
        
        # Integrate current phase
        problem = _ode_problem(mps.phases[k], (t_current, t_end), z_current, v)
        raw_sol = _solve(problem, get_alg(mps), get_options(mps))
        
        # Extract final state
        z_current = raw_sol[end]
        
        # Apply jump if not at final time
        if t_end < tf && k < length(mps.phases)
            Δx, Δp = mps.transitions[k].jump
            z_current = apply_jump(z_current, Δx, Δp)
        end
        
        t_current = t_end
        
        # Exit if we've reached tf
        if t_current >= tf
            break
        end
    end
    
    return split_state_costate(z_current)  # (xf, pf)
end

Complete Solution Logic

For complete solution sol = f((t0, tf), x0, p0, v):

function evaluate_solution(mps::MultiPhaseSystem, tspan, z0, v)
    t0, tf = tspan
    phase_idx = find_phase(mps, t0)
    
    # Solve each phase
    solutions = []
    z_current = z0
    t_current = t0
    
    for k in phase_idx:length(mps.phases)
        t_end = k < length(mps.phases) ? mps.transitions[k].time : tf
        t_end = min(t_end, tf)
        
        # Solve phase
        problem = _ode_problem(mps.phases[k], (t_current, t_end), z_current, v)
        phase_sol = _solve(problem, get_alg(mps), get_options(mps))
        push!(solutions, phase_sol)
        
        # Prepare for next phase
        if t_end < tf && k < length(mps.phases)
            z_current = phase_sol[end]
            Δx, Δp = mps.transitions[k].jump
            z_current = apply_jump(z_current, Δx, Δp)
        end
        
        t_current = t_end
        
        if t_current >= tf
            break
        end
    end
    
    # Merge solutions (delegated to extension)
    return _merge_solutions(solutions, mps)
end

Solution Merging (Extension)

The merging logic is delegated to the SciMLBase extension, following the pattern from DiffEqParamEstim.jl:

# ext/CTFlowsSciMLBaseExt.jl
function CTFlows._merge_solutions(solutions::Vector, mps::CTFlows.MultiPhaseSystem)
    # Concatenate time and state vectors, removing duplicates at transitions
    u = [uc for (k, sol) in enumerate(solutions) for uc in sol.u[1:end-1]]
    t = [tc for (k, sol) in enumerate(solutions) for tc in sol.t[1:end-1]]
    
    # Add final point
    push!(u, solutions[end].u[end])
    push!(t, solutions[end].t[end])
    
    # Build merged solution using SciMLBase
    # Note: We create a custom solution type that holds the phase information
    merged = CTFlows.MergedSolution(u, t, solutions)
    
    # Use SciMLBase to build a proper ODESolution
    return SciMLBase.build_solution(
        solutions[1].prob,  # Use first problem as template
        solutions[1].alg,
        merged.t,
        merged.u,
        retcode = :Success
    )
end

---

3.4 Summary of Internal Types

Role	Type/Function	Purpose
Flow	Flow{S<:AbstractSystem}	User-facing callable, stores system + options
System	AbstractSystem{T}	Encapsulates dynamics, parameterized by type T (e.g., :state_feedback)
Multi-Phase	MultiPhaseSystem{S}	Manages concatenated flows with phase transitions
Problem Builder	_ode_problem(system, ...)	Stub in src/, overridden in extension
Solver	_solve(problem, alg, options)	Stub in src/, overridden in extension
Solution Merger	_merge_solutions(solutions, mps)	Stub in src/, overridden in extension
Post-processor	build_solution(raw_sol, system)	Transforms raw solution to typed result

---

3.5 Pseudo-code Call Graphs

1. Simple Flow Creation

# User API
f = Flow(ocp, u, :state_feedback, alg=Tsit5(), abstol=1e-8)

# Internal Factory (src/)
system = build_system(ocp, u, :state_feedback)   # Compose RHS
options = (alg=Tsit5(), abstol=1e-8, ...)
return Flow(system, options)

2. Point Evaluation

# User API
xf, pf = f(t0, x0, p0, tf, v; reltol=1e-6)

# Internal (src/)
z0 = (x0, p0)
merged_options = merge(get_options(f), (reltol=1e-6,))

# Stub (overridden in extension)
problem = _ode_problem(get_system(f), (t0, tf), z0, v)
raw_sol = _solve(problem, merged_options.alg, merged_options)

return get_final_state(raw_sol, get_system(f))  # (xf, pf)

3. Complete Solution

# User API
sol = f((t0, tf), x0, p0, v)

# Internal (src/)
z0 = (x0, p0)
problem = _ode_problem(get_system(f), (t0, tf), z0, v)
raw_sol = _solve(problem, get_options(f).alg, get_options(f))

return build_solution(raw_sol, get_system(f))  # OptimalControlSolution

4. Concatenation (Point)

# User API
f1 = Flow(ocp, u1, :state_feedback)
f2 = Flow(ocp, u2, :state_feedback)
f_concat = f1 * (t1, Δx, Δp, f2)

# Internal (src/)
# Type check: ensure get_flow_type(f1) == get_flow_type(f2)
mps = MultiPhaseSystem(
    phases = [get_system(f1), get_system(f2)],
    transitions = [PhaseTransition(t1, (Δx, Δp))]
)
return Flow(mps, get_options(f1))

# At call time: xf, pf = f_concat(t0, x0, p0, tf, v)
# Calls evaluate_point(mps, t0, (x0, p0), tf, v)

5. Concatenation (Solution)

# User API
sol = f_concat((t0, tf), x0, p0, v)

# Internal (src/)
# Calls evaluate_solution(mps, (t0, tf), (x0, p0), v)
# Which:
#   1. Solves each phase
#   2. Applies jumps between phases
#   3. Calls _merge_solutions(solutions, mps) [extension]

---

4. Design Decisions (V3)

1. No Wrappers: HamiltonianFeedback, StateFeedback, etc., are removed from the public API.
2. No Descriptions: Symbols like :hamiltonian are not used for dispatch in the factory; the type of the first argument and positional control argument are sufficient.
3. Positional Control Argument: Control specification (:state_feedback, :dynamic_feedback, etc.) is a positional argument, not a named argument. Named arguments are reserved for solver options.
4. Duck-typed Initial Conditions: Internal rhs! must be robust to the shape of the input state.
5. Lightweight Options: Simple NamedTuple merging, no OCPTool overhead.
6. Extension-Based Backend: SciMLBase extension provides ODEProblem, solve, and solution merging. Stubs in src/ provide helpful error messages.
7. Type-Safe Concatenation: Flows can only be concatenated if they are of the same strict type. For OCP flows, this means the same control mode (tracked via type parameter AbstractSystem{T}).
8. Multi-Phase Architecture: Concatenated flows are represented as MultiPhaseSystem with explicit phase transitions and jumps. Solution merging follows the DiffEqParamEstim.jl pattern.
9. No Hamiltonian Flows: V3 focuses on OCP and ODEFunction flows only. Hamiltonian flows are removed to simplify the architecture.

---

5. Advantages of V3 Architecture

1. Minimal Dependencies: Core package has no heavy dependencies. SciMLBase is loaded only when needed.
2. Clear Error Messages: Stubs provide actionable guidance when extension is not loaded.
3. Efficient Concatenation: Multi-phase system avoids creating intermediate Flow objects.
4. SciML Compatibility: Solution merging reuses proven patterns from the SciML ecosystem.
5. Type Safety: Compile-time guarantees prevent invalid flow concatenations.
6. Extensible: Easy to add new system types or backends in the future.
7. Performance: Lightweight core, efficient phase transitions, minimal allocations.

---

6. Migration from V2

Key changes from V2:

1. Removed: Hamiltonian flows (use OCP flows instead)
2. Changed: Control specification is now positional (Flow(ocp, u, :state_feedback) instead of control=:state_feedback)
3. Changed: Concatenation creates MultiPhaseSystem instead of nested Flow objects
4. Changed: Extension triggered by SciMLBase instead of OrdinaryDiffEq
5. Added: Explicit stub functions with helpful error messages
6. Added: Solution merging logic following DiffEqParamEstim.jl pattern

---

7. Future Considerations

1. GPU Support: Multi-phase architecture is compatible with GPU arrays
2. Parallel Shooting: MultiPhaseSystem can be extended to support parallel phase integration
3. Adaptive Phase Refinement: Automatic subdivision of phases based on error estimates
4. Custom Backends: While V3 focuses on SciML, the stub pattern allows for alternative backends
5. Sensitivity Analysis: Phase structure enables efficient adjoint computation

[jbcaillau]

@ocots priority on

• flows from ocp
• performance

[ocots]

CTFlows.jl V3 Refactoring Planning

Issue: #143 - [Dev] Refactoring
Date: 2025-12-18
Status: Planning Complete ✅

TL;DR

Refactor CTFlows.jl to V3 architecture: remove type wrappers from differential geometry tools, simplify Flow API to OCP and ODEFunction only, implement SciMLBase extension with helpful stubs, and introduce multi-phase concatenation with solution merging. This will result in a simpler, more performant, and more maintainable codebase aligned with the V3 architecture report.

---

1. Overview

Goal

Implement the V3 architecture for CTFlows.jl as documented in reports/architecture_use_cases_v3.md. This involves:

• Removing type wrappers (Hamiltonian, VectorField, HamiltonianLift) from differential geometry tools
• Simplifying Flow API to support only OCP and ODEFunction flows
• Implementing SciMLBase extension pattern with helpful error messages
• Introducing multi-phase concatenation with proper solution merging

Key Features

• Wrapper-Free Differential Geometry: All tools work with pure Function objects
• Simplified Flow API: Only OCP and ODEFunction flows (no standalone Hamiltonian flows)
• SciMLBase Extension: Lightweight core with extension-based backend integration
• Multi-Phase Concatenation: Proper handling of concatenated flows with jumps and solution merging
• Positional Control Argument: Control specification as positional argument (not named)

References

• V3 Architecture Report
• Goddard Tutorial (reference implementation)
• DiffEqParamEstim.jl (solution merging pattern)

---

2. User Stories

ID	Description	Status
US-1	As a user, I want to use differential geometry tools without managing type wrappers (Hamiltonian, VectorField, HamiltonianLift)	⏳
US-2	As a user, I want to create flows from OCPs using positional control arguments (:state_feedback, :dynamic_feedback, :open_loop) instead of type wrappers (ControlLaw, FeedbackControl)	⏳
US-3	As a user, I want to create flows from ODEFunctions (OOP and IP)	⏳
US-4	As a user, I want to concatenate flows with automatic phase management	⏳
US-5	As a user, I want clear error messages when SciMLBase is not loaded	⏳
US-6	As a user, I want automatic variable inference from OCP properties (already implemented ✅)	⏳

---

3. Technical Decisions

Decision	Choice
Type Wrappers	Remove all wrappers (Hamiltonian, VectorField, HamiltonianLift)
Control Argument	Positional (:state_feedback, :dynamic_feedback, etc.) not named
Backend Extension	SciMLBase extension with stubs in src/
Concatenation	MultiPhaseSystem with explicit phase transitions
Solution Merging	Follow DiffEqParamEstim.jl pattern in extension
Hamiltonian Flows	Remove standalone Hamiltonian flows (use OCP instead)
Test Infrastructure	Standard Pkg.test (no CTBase.run_tests)
Coverage	Manual coverage via julia --project=. --code-coverage=user

---

4. Tasks

Phase 1: Differential Geometry Refactoring

Task	Description
T1.1	Remove VectorField, Hamiltonian, HamiltonianLift type definitions from src/types.jl
T1.2	Refactor Lift(f) to return pure Function instead of HamiltonianLift
T1.3	Remove ⋅ operator (cannot detect autonomy without wrappers)
T1.4	Keep only Lie(X, f; autonomous, variable) with explicit arguments
T1.5	Refactor Lie bracket: rename to ad(X, Y; autonomous, variable) to avoid ambiguity with Lie derivative
T1.6	Refactor Poisson(H, G; autonomous, variable) to work with pure functions
T1.7	Update @Lie macro: remove Hamiltonian(...) wrapper creation, expand @Lie [X, Y] to ad(X, Y) instead of Lie(X, Y)
T1.8	Update all differential geometry tests in test/test_differential_geometry.jl

Phase 2: Flow API Simplification

Task	Description
T2.1	Remove standalone Hamiltonian flow constructors (ext/hamiltonian.jl)
T2.2	Remove VectorField flow constructors (ext/vector_field.jl)
T2.3	Remove ControlLaw, FeedbackControl, StateConstraint, MixedConstraint, Multiplier wrappers from src/types.jl
T2.4	Create new OCP flow constructors with optional control_type (default :dynamic_feedback): Flow(ocp, u, control_type::Symbol = :dynamic_feedback; alg, ...) Flow(ocp, u, control_type::Symbol, g, μ; alg, ...) (with constraints)
T2.5	Implement control type dispatch logic in ext/optimal_control_problem.jl
T2.6	Update ODEFunction flow constructor (OOP and IP support)
T2.7	Update all flow tests (test/test_flow_*.jl, test/test_optimal_control_problem.jl)
T2.8	Remove test/test_flow_hamiltonian.jl, test/test_flow_vector_field.jl, test/test_flow_hamiltonian_vector_field.jl

Phase 3: Extension and Stub Implementation

Task	Description
T3.1	Create src/stubs.jl with _ode_problem and _solve stubs
T3.2	Implement helpful error messages in stubs
T3.3	Rename ext/CTFlowsODE.jl to ext/CTFlowsSciMLBaseExt.jl
T3.4	Update extension to override _ode_problem and _solve
T3.5	Update Project.toml with SciMLBase extension
T3.6	Test extension loading and error messages

Phase 4: Multi-Phase Concatenation

Task	Description
T4.1	Create src/concatenation.jl with AbstractMultiPhaseSystem{S}, MultiPhaseSystem{S}, and PhaseTransition types
T4.2	Implement evaluate_point(mps::AbstractMultiPhaseSystem, ...) in src/ (integrates each phase sequentially, applies jumps)
T4.3	Implement evaluate_solution(mps::AbstractMultiPhaseSystem, ...) in src/ (stores all phase solutions for merging)
T4.4	Create _merge_solutions(solutions::Vector, mps::AbstractMultiPhaseSystem) stub in src/ with helpful error message
T4.5	Implement _merge_solutions in ext/CTFlowsSciMLBaseExt.jl using batch merging strategy (concatenate u/t vectors, use SciMLBase.build_solution)
T4.6	Update concatenation operator * to create MultiPhaseSystem instead of concatenating RHS
T4.7	Update concatenation tests in test/test_concatenation.jl to test multi-phase evaluation and solution merging

Phase 5: Documentation and Migration

Task	Description
T5.1	Update README.md with V3 examples
T5.2	Create migration guide (V2 → V3)
T5.3	Update docstrings for all public API functions
T5.4	Update Goddard tutorial to use V3 API
T5.5	Add V3 architecture report to documentation

---

5. Testing Guidelines

Test file structure

# test/test_<feature>.jl

function test_<feature>()
    # ========================================================
    # Unit tests – helper logic (no side effects)
    # ========================================================
    @testset "<feature> helpers" begin
        # Test individual helper functions
    end
    
    # ========================================================
    # Integration tests – actual functionality
    # ========================================================
    @testset "<feature> integration" begin
        # Test full feature with real OCP/ODEFunction
    end
end

Key test files to update

• test/test_differential_geometry.jl (Phase 1)
• test/test_flow_hamiltonian.jl (Phase 2 - remove)
• test/test_flow_vector_field.jl (Phase 2 - remove)
• test/test_flow_hamiltonian_vector_field.jl (Phase 2 - remove)
• test/test_optimal_control_problem.jl (Phase 2)
• test/test_flow_function.jl (Phase 2)
• test/test_concatenation.jl (Phase 4)

---

6. Test Commands

# Run all tests
julia --project=. -e 'using Pkg; Pkg.test("CTFlows");'

# Run specific test file
julia --project=. test/test_differential_geometry.jl

# Run with coverage
julia --project=. --code-coverage=user -e 'using Pkg; Pkg.test("CTFlows");'

---

7. Coverage Testing

Coverage command

julia --project=. --code-coverage=user -e 'using Pkg; Pkg.test("CTFlows");'

Target

≥ 90% coverage for refactored code.

Iteration process

1. Run coverage command
2. Check uncovered lines: find src -name "*.jl.*.cov"
3. Add tests for uncovered code paths
4. Repeat until ≥ 90%

---

8. GitHub Workflow

Structure

Issue #143 ([Dev] Refactoring) ← Epic
  ├── PR "Phase 1: Differential Geometry Refactoring" → linked to #143
  ├── PR "Phase 2: Flow API Simplification" → linked to #143
  ├── PR "Phase 3: Extension and Stub Implementation" → linked to #143
  ├── PR "Phase 4: Multi-Phase Concatenation" → linked to #143
  └── PR "Phase 5: Documentation and Migration" → closes #143

Checklist for Issue #143

• Phase 1: Differential Geometry Refactoring
• Phase 2: Flow API Simplification
• Phase 3: Extension and Stub Implementation
• Phase 4: Multi-Phase Concatenation
• Phase 5: Documentation and Migration

---

9. MVP (Minimum Viable Product)

MVP = Phase 1 + Phase 2 + Phase 3

The MVP includes:

• Wrapper-free differential geometry tools (all existing functionality preserved)
• Simplified Flow API (OCP and ODEFunction only)
• SciMLBase extension with helpful error messages
• All existing tests passing with updated API

Not in MVP (deferred to later phases):

• Multi-phase concatenation (Phase 4)
• Complete documentation update (Phase 5)

---

10. Breaking Changes

API Changes

Old API (V2 - Current)	New API (V3 - Target)	Migration
u = ControlLaw((x,p) -> p[2]); Flow(ocp, u)	Flow(ocp, (x,p) -> p[2]) or Flow(ocp, (x,p) -> p[2], :dynamic_feedback)	Remove ControlLaw wrapper; control_type optional (default :dynamic_feedback)
u = FeedbackControl(x -> -x); Flow(ocp, u)	Flow(ocp, x -> -x, :state_feedback)	Remove FeedbackControl wrapper, specify :state_feedback
Flow(ocp, u; autonomous=false)	Flow(ocp, u, :dynamic_feedback, autonomous=false)	Add optional control_type before named args
Flow(ocp, u, g, μ) (with wrappers)	Flow(ocp, u, :dynamic_feedback, g, μ)	Remove wrappers, add control_type after u
H = Hamiltonian(f)	H = f (pure function)	Remove wrapper
X = VectorField(f)	X = f (pure function)	Remove wrapper
HL = HamiltonianLift(X)	HL = Lift(X) (returns Function)	Use Lift function
Flow(H::Hamiltonian)	Flow(ocp, ...) (use OCP instead)	Convert to OCP flow
Flow(X::VectorField)	Flow(ode_function) (use ODEFunction)	Convert to ODEFunction flow

Deprecation Strategy

1. Phase 1-3: Implement new API alongside old (if feasible)
2. Phase 4: Remove old API completely
3. Phase 5: Document migration path

---

11. Risk Assessment

Risk	Impact	Mitigation
Breaking existing user code	High	Provide detailed migration guide
Test coverage gaps	Medium	Aim for ≥ 90% coverage
Extension loading issues	Medium	Comprehensive stub error messages
Solution merging bugs	High	Follow proven DiffEqParamEstim pattern
Performance regression	Low	Benchmark critical paths

---

12. Success Criteria

• All tests pass with V3 API
• ≥ 90% code coverage
• Goddard tutorial works with V3 API
• No type wrappers in differential geometry
• SciMLBase extension loads correctly
• Clear error messages when extension not loaded
• Migration guide complete
• Documentation updated