Skip to content

Conversation

@jbcaillau
Copy link
Member

@jbcaillau jbcaillau commented Dec 8, 2025

start with going to x[i1, i2 in 1:1, k in 1:1, j] for M(n, 1, R)

  • 0:grid_size to 1:grid_size+1 to avoid discrepancy between x and x_m
  • [ ] implement operations on x_m instead of x
  • at this step, add idoneus subs to replace x(t) or x[rg](g) by x[rg, j] (should work)
  • same update can be made for u to allow u(t) -> u[:, j], etc., still with (real) vectors
  • vectorise dynamics for exa
  • vectorise constraints for exa
  •  at this step, re-write some qc tests (with tweaks / macros for matrix / complex values)
  • go for tensors (matrix + complex valued)
  • for that, use the fact that Julia vectors are indeed n x 1 matrices (v[i] == v[i, 1] for a Vector), consistently with R^n = M(n, 1, R)

NB:

  • needed overloads on exa : *, adjoint, dot... (some work, e.g. sum)
  • use @def_exa to test
  • need some :fun stuff to plot, etc. if going to tensors

@jbcaillau jbcaillau linked an issue Dec 8, 2025 that may be closed by this pull request
@github-actions
Copy link
Contributor

github-actions bot commented Dec 8, 2025

Breakage test results
Date: 2025-12-13 17:36:38

Name Latest Stable
OptimalControl compat: v0.7.2 compat: v0.7.2

@jbcaillau jbcaillau changed the title first commit pr-adding-tensors Dec 8, 2025
jbcaillau and others added 7 commits December 8, 2025 18:23
Change all 0:grid_size indexing to 1:grid_size+1 in the exa backend:
- State and control variable creation ranges
- Dynamics constraint loops (0:grid_size-1 → 1:grid_size)
- Lagrange cost integration loops with proper scheme handling
- Path constraint loops (0:grid_size → 1:grid_size+1)
- Initial/final condition indexing (0/grid_size → 1/grid_size+1)
- Mayer cost boundary state indexing

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>
Replace :grid_size+1 with :(grid_size+1) for proper Julia expression
quoting in final state index arguments:
- Boundary constraints (subs2 call)
- Final constraints (subs3 call)
- Mayer cost (subs2 call)

Keeps range expressions like 1:grid_size+1 unchanged.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>
This migration enables standard matrix operations (range slicing, arithmetic,
broadcasting) on state variables in the exa backend by using p.x_m (a matrix
of ExaModels.Var) instead of p.x (ExaModels.Variable).

Changes:
- Add x_m field to ParsingInfo struct
- Initialize p.x_m in p_state! and p_state_exa! functions
- Replace p.x with p.x_m in 13 locations across 4 functions:
  * p_constraint_exa!: boundary and path constraints (5 replacements)
  * p_dynamics_coord_exa!: state dynamics discretization (4 replacements)
  * p_lagrange_exa!: Lagrange cost function (2 replacements)
  * p_mayer_exa!: Mayer cost function (2 replacements)
- Remove 3 useless substitution lines in p_constraint_exa!
- Enable direct arithmetic operations: p.x_m[i, j+1] - p.x_m[i, j]

The p.x_m matrix is defined as [x[i, j] for i ∈ 1:n, j ∈ 1:grid_size+1],
providing more efficient access patterns for discretized optimal control.

Related to issue #181

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Ignore Claude context files used for development documentation that
should not be tracked in version control.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
@jbcaillau
Copy link
Member Author

jbcaillau commented Dec 9, 2025

issue: x_m[i, j] ... for j in ... with j that is expected by ExaModels to generate a generator (inside dynamics, constraints...), not an array, as generators (lazy eval, no alloc) are used for SIMD modelling.

fix: don't use x_m but x[i, j] with j generator and i not a generator but a standard index to define a comprehension (then a vector) as below.

so:

  • come back to p.x everywhere
  • add subs to replace x(t), x[rg](t) by [x[i, j] for i in 1:n / rg] for dynamics, constraints, and costs
julia> g(x) = [x[1] * x[2] + x[3], x[3] * x[4]]
g (generic function with 1 method)

julia> x
Variable

  x  R^{4 × 10}


julia> constraint(c, g([x[i, j] for i in 1:4])[1] for j in 1:10)
Constraint

  s.t. (...)
       g♭  [g(x,θ,p)]_{p  P}  g♯

  where |P| = 10


julia> for k in 1:2 
       constraint(c, g([x[i, j] for i in 1:4])[k] for j in 1:10)
       end

@jbcaillau
Copy link
Member Author

non scalar dynamics, (nonlinear) constraints:

  • use a loop as above
  • if the constraint is scalar, use the standard ct convention and treat is as scalar (not length one vector)

jbcaillau and others added 3 commits December 12, 2025 10:25
This commit addresses the p.x_m indexing bug and standardizes on 0-based
grid indexing for the ExaModels backend.

Changes:
1. Removed p.x_m matrix wrapper from ParsingInfo struct
2. Reverted all p.x_m references back to p.x (ExaModels.Variable)
   - p.x properly handles symbolic indexing in ExaModels generators
   - p.x_m (plain Matrix) was causing ArgumentError with symbolic indices

3. Migrated from 1-based to 0-based grid indexing:
   - State/control variables: 0:grid_size (was 1:grid_size+1)
   - Initial conditions: index 0 (was 1)
   - Final conditions: index grid_size (was grid_size+1)
   - Dynamics constraints: 0:grid_size-1 (was 1:grid_size)
   - Path constraints: 0:grid_size (was 1:grid_size+1)
   - Lagrange cost integration:
     * Euler forward: 0:grid_size-1 (was 1:grid_size)
     * Euler implicit: 1:grid_size (was 2:grid_size+1)
     * Midpoint: 0:grid_size-1 (was 1:grid_size)
     * Trapeze: endpoints (0, grid_size), interior 1:grid_size-1

4. Documentation improvements in utils.jl for subs3 and subs4

Rationale:
The p.x_m approach attempted to enable matrix operations but failed because
generator variables in ExaModels become symbolic types (ParSource, Node2),
which cannot index regular Julia arrays. Using p.x (ExaModels.Variable)
throughout ensures proper symbolic indexing support.

0-based indexing provides cleaner semantics for discretization schemes,
with grid points at times t0 + j*dt for j in 0:grid_size.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
@jbcaillau
Copy link
Member Author

jbcaillau commented Dec 12, 2025

  • CRITICAL: review all $p.foo$ vs. $(p.foo)$ expressions in onepass.jl
  • as a rule: $(p.foo) when the value is a symbol that needs to be interpretated in the (escaped) call context

FTR

julia> n = 6
6

julia> f.x = :n
:n

julia> e = :($(f.x) == 6)
:(n == 6)

julia> eval(e)
true

julia> e = :($f.x == 6)
:((Foo(:n)).x == 6)

julia> eval(e)
false

julia> e = :($(f.x) == 6)
:(n == 6)

julia> f.x = 6
6

julia> e = :($(f.x) == 6)
:(6 == 6)

… (min/max) as these two things are known statically (= ater parsing, without evaluating), plus the fact p.criterion is a symbol (:min...)
Relocated these utility functions to utils.jl for better organization:
- is_range(x): Predicate to test if x is a range (AbstractRange or i:j expr)
- as_range(x): Normalizes x to a range or single-element array

Also added clarifying comments at all as_range call sites explaining
the "case rg = i (vs i:j or i:p:j)" normalization.

These functions are general-purpose utilities used for constraint
index handling and belong with other utility functions rather than
in the parsing logic.
Major enhancements:

1. Enhanced subs2 function (utils.jl):
   - Added range indexing support: x[1:3] → [y[k, j] for k ∈ 1:3]
   - Preserves scalar indexing behavior: x[i] → y[i, j]
   - Added optional 'k' parameter for predictable symbol generation in tests
   - New pattern: :($xx[$rg]) when is_range(rg) generates comprehension
   - Existing pattern: :($xx[$i]) continues to work for scalars

2. Reorganized utility functions (onepass.jl → utils.jl):
   - Moved __symgen() to utils.jl (needed by subs2 for range support)
   - Consolidated is_range() and as_range() at top of utils.jl
   - Better organization: all utility functions now in one place

3. Comprehensive test suite (test_utils.jl):
   - 18 tests for subs2 (up from 2 original tests)
   - Tests organized into 6 testsets:
     * Scalar indexing (backward compatibility - 4 tests)
     * Range indexing (new functionality - 5 tests)
     * Mixed scalar/range (1 test)
     * Nested/complex expressions (3 tests)
     * Edge cases (4 tests)
     * Backward compatibility verification (1 test)
   - All tests use explicit k parameter for exact equality assertions

4. Test configuration (runtests.jl):
   - Default to running only utils and utils_bis tests
   - Other tests disabled to speed up development iteration
   - Can still run all tests with: julia test/runtests.jl all

Implementation details:
- Range detection uses is_range() predicate (checks AbstractRange or i:j expr)
- Comprehension variable k generated via __symgen() or passed explicitly
- Backward compatible: existing code using scalar indices unchanged
- Forward compatible: supports symbolic ranges (1:n, 2:2:grid_size, etc.)

This enhancement enables tensor/matrix operations in the ExaModels backend
by allowing range-based substitutions in symbolic expressions.
Re-enabled all tests (aqua, prefix, onepass_fun, onepass_exa) that were
temporarily disabled during subs2 development and testing.
@jbcaillau
Copy link
Member Author

jbcaillau commented Dec 13, 2025

for vectorised access to state in exa:

  • treat x(t) by improving subs2
  • update tests (and docstrings), both for subs2...
  • ... and dynamics (sum on whole / ranges for state / control, call to user defined functions...)
  • same for subs5 (and rename it subs2m)
  • remove subs4 from utils tests
  • same for constraints and costs

jbcaillau and others added 3 commits December 13, 2025 12:55
This commit completes the migration to the enhanced subs2 function that now
handles three patterns: scalar indexing, range indexing, and bare symbols.

Changes:
- src/onepass.jl: Updated all 14 subs2 calls to include the new dim parameter
  - State substitutions use p.dim_x
  - Control substitutions use p.dim_u
  - Locations: p_constraint_exa! (lines 700-801), p_dynamics_coord_exa! (lines 900-907),
    p_lagrange_exa! (lines 970-974), p_mayer_exa! (lines 1033-1034)

- test/test_utils.jl: Updated all subs2 test calls to include dim parameter
  - Modified test expectations for bare symbol behavior (now expands to comprehensions)
  - Added explicit k parameter where needed for predictable test results
  - Updated backward compatibility tests to reflect new bare symbol pattern

- src/utils.jl: Enhanced subs2 docstring
  - Documents all three substitution patterns
  - Explains the dim parameter for bare symbol expansion
  - Provides comprehensive examples for scalar, range, and bare symbol usage

The new subs2 pattern `x` (bare symbol) → `[y[k, j] for k ∈ 1:dim]` enables
tensor operations where entire state/control vectors are referenced without
explicit indexing.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
This commit reverts the dim parameter addition to subs2, as bare symbol
expansion (x → [y[k, j] for k ∈ 1:dim]) cannot be handled at the subs2 level.

Changes:
- src/onepass.jl: Removed dim parameter from all 14 subs2 calls
  - Lines 700-701: Boundary constraints (p_constraint_exa!)
  - Lines 800-801: Path constraints (p_constraint_exa!)
  - Lines 900-907: Dynamics discretization (p_dynamics_coord_exa!)
  - Lines 970-974: Lagrange cost (p_lagrange_exa!)
  - Lines 1033-1034: Mayer cost (p_mayer_exa!)

- src/utils.jl: Updated subs2 docstring
  - Removed bare symbol pattern from documentation
  - Changed from 3 patterns to 2 patterns (scalar and range only)
  - Updated examples to show bare symbols are NOT substituted
  - Restored original signature: subs2(e, x, y, j; k = __symgen(:k))

- test/test_utils.jl: Updated all subs2 tests
  - Removed dim parameter from all 18 test calls
  - Updated Test 2 and Test 19 expectations (bare symbols NOT substituted)
  - Verified all range indexing tests still pass without dim parameter

Rationale: Bare symbol expansion requires dimension information that is not
available at the expression substitution level. This functionality will need
to be implemented at a higher level in the parsing pipeline.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Implements the pattern from p_constraint_exa! boundary case across all functions
that use subs2, enabling bare symbol handling (e.g., x(t) → [x[k, j] for k ∈ 1:dim_x]).

Changes to src/onepass.jl:
1. p_mayer_exa! (lines 1030-1044): +3 lines
   - Added k = __symgen(:k)
   - Added subs for x0 and xf bare symbols

2. p_constraint_exa! path constraints (lines 797-815): +3 lines
   - Added k = __symgen(:k)
   - Added subs for xt and ut bare symbols

3. p_lagrange_exa! (lines 968-985): +4 lines
   - Added k = __symgen(:k)
   - Added subs for xt and ut bare symbols in ej1 and ej12

4. p_dynamics_coord_exa! (lines 900-920): +6 lines
   - Added k = __symgen(:k)
   - Added subs for xt and ut bare symbols in ej1, ej2, and ej12

Total: 16 lines added across 4 functions

Pattern applied:
- subs2 handles indexed cases: x[i] → y[i, j] and x[1:3] → [y[k, j] for k ∈ 1:3]
- subs handles bare symbols: x → [y[k, j] for k ∈ 1:dim]
- Same symbol k used for both state and control (sequential application)

This completes the implementation of tensor support for bare symbols in the
ExaModels backend.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
@jbcaillau
Copy link
Member Author

for vectorised access to state in exa:

  • treat x(t) by combiningsubs2 and subs
  • update tests (and docstrings), both for subs2 and dynamics
  • same for subs5 (and rename it subs2m)
  • remove subs4 from utils tests
  • same for constraints and costs

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[Dev] Adding tensors

2 participants