test: deterministic seeds + multi-seed robustness sweeps for flaky tests

README.md

@@ -100,8 +100,38 @@ solve!(prob; max_iter=100)
 - **Modular objectives**: Combine cost terms with `+` and `*` (regularization, minimum time, terminal cost, etc.)
 - **Constraint support**: Bounds, equality, nonlinear, symmetry, and L1 slack constraints
 - **Automatic differentiation**: Sparse Jacobians and Hessians via ForwardDiff
+- **Sparse formulations**: Exploits problem structure for efficiency
 - **Solver callbacks**: Monitor and control the optimization process
 
+## Testing
+
+```bash
+julia --project=. test/runtests.jl
+```
+
+`runtests.jl` runs every `@testitem` in `src/`, `ext/`, and `test/`. Tests in
+`benchmark/` are skipped — that subdirectory ships its own `Project.toml`
+(extra deps like `HarmoniqsBenchmarks`) and has a dedicated workflow.
+
+### Stochastic / numerical primitives — two-layer testing
+
+A single seeded `MersenneTwister` is reproducible on one Julia version but
+small downstream numerics can drift across the CI matrix (1.10 / 1.11 / 1.12).
+For tests that touch non-deterministic surfaces (solver convergence from
+random init, finite-difference derivative comparisons) we pair each test:
+
+1. **Deterministic baseline**: a single seeded trajectory + multiplier.
+   A failure is a real regression on a specific (Julia version, seed) pair.
+2. **Robustness sweep**: K=20 independent seeds; passes if a fraction of
+   seeds (per-test threshold, chosen with buffer above the observed baseline
+   rate — typically 0.80, lower for inherently noisy checks like norm-based
+   finite-diff) land within tolerance. Detects regressions that drop the
+   true pass rate well below the threshold with very high probability
+   (binomial), while staying insensitive to lucky/unlucky single draws.
+
+The sweeps are cheap enough (a handful of seconds in aggregate) to run on
+every PR.
+
 ## Contributing
 
 ### Building Documentation

src/constraints/_constraints.jl

@@ -24,6 +24,7 @@ using FiniteDiff
 using SparseArrays
 using TestItemRunner
 using LinearAlgebra
+using Random
 using Test
 
 # Import and extend the common interface
@@ -121,7 +122,8 @@ function get_full_hessian end
         show_hessian_diff=false,
         test_equality=true,
         atol=1e-5,
-        rtol=1e-5
+        rtol=1e-5,
+        rng=Random.default_rng(),
     )
 
 Test that constraint Jacobian and Hessian match finite difference approximations.
@@ -136,9 +138,16 @@ Test that constraint Jacobian and Hessian match finite difference approximations
 - `test_equality=true`: Test element-wise equality (vs norm-based test)
 - `atol=1e-5`: Absolute tolerance
 - `rtol=1e-5`: Relative tolerance
+- `rng`: RNG used to draw the Lagrange multiplier `μ` for the Hessian check.
+  Pass a seeded `AbstractRNG` (e.g. `MersenneTwister(seed)`) for reproducible
+  test results across runs and Julia versions.
+- `assert=true`: When `true`, fails the surrounding `@testset` via `@test` on
+  mismatch. Set `false` to inspect `jacobian_pass`/`hessian_pass` in the
+  returned NamedTuple without registering a test outcome — useful for
+  multi-seed robustness sweeps.
 
 # Returns
-Tuple of (∂g, ∂g_finite_diff, μ∂²g, μ∂²g_finite_diff) for inspection
+NamedTuple `(; ∂g, ∂g_finite_diff, μ∂²g, μ∂²g_finite_diff, jacobian_pass, hessian_pass)`
 
 # Example
 ```julia
@@ -155,6 +164,8 @@ function test_constraint(
     test_equality = true,
     atol = 1e-5,
     rtol = 1e-5,
+    rng::AbstractRNG = Random.default_rng(),
+    assert::Bool = true,
 )
 
     # Function to evaluate constraint via evaluate!
@@ -192,18 +203,19 @@ function test_constraint(
     end
 
     # Test Jacobian equality
-    if test_equality
-        @test all(isapprox.(∂g, ∂g_finite_diff, atol = atol, rtol = rtol))
+    jac_pass = if test_equality
+        all(isapprox.(∂g, ∂g_finite_diff, atol = atol, rtol = rtol))
     else
         if atol > 0.0
-            @test norm(∂g - ∂g_finite_diff) < atol
+            norm(∂g - ∂g_finite_diff) < atol
         else
-            @test norm(∂g - ∂g_finite_diff) / norm(∂g_finite_diff) < rtol
+            norm(∂g - ∂g_finite_diff) / norm(∂g_finite_diff) < rtol
         end
     end
+    assert && @test jac_pass
 
     # Test Hessian
-    μ = rand(constraint.dim)
+    μ = rand(rng, constraint.dim)
 
     μ∂²g = CommonInterface.eval_hessian_of_lagrangian(constraint, traj, μ)
 
@@ -225,17 +237,25 @@ function test_constraint(
     end
 
     # Test Hessian equality (only upper triangle since Hessian is symmetric)
-    if test_equality
-        @test all(isapprox.(triu(μ∂²g), triu(μ∂²g_finite_diff), atol = atol))
+    hess_pass = if test_equality
+        all(isapprox.(triu(μ∂²g), triu(μ∂²g_finite_diff), atol = atol))
     else
         if atol > 0.0
-            @test norm(μ∂²g - μ∂²g_finite_diff) < atol
+            norm(μ∂²g - μ∂²g_finite_diff) < atol
         else
-            @test norm(μ∂²g - μ∂²g_finite_diff) / norm(μ∂²g_finite_diff) < rtol
+            norm(μ∂²g - μ∂²g_finite_diff) / norm(μ∂²g_finite_diff) < rtol
         end
     end
-
-    return ∂g, ∂g_finite_diff, μ∂²g, μ∂²g_finite_diff
+    assert && @test hess_pass
+
+    return (;
+        ∂g,
+        ∂g_finite_diff,
+        μ∂²g,
+        μ∂²g_finite_diff,
+        jacobian_pass = jac_pass,
+        hessian_pass = hess_pass,
+    )
 end
 
 export test_constraint

src/constraints/linear/time_consistency_constraint.jl

@@ -86,17 +86,23 @@ end
 @testitem "TimeConsistencyConstraint with free time optimization" begin
     include("../../../test/test_utils.jl")
     using NamedTrajectories
+    using Random
+
+    # Deterministic seed: same trajectory + multipliers on every run. A failure
+    # here is a real regression, not RNG drift. Robustness across seeds is
+    # covered by the `:robustness` testitem below.
+    rng = MersenneTwister(0)
 
     # Create trajectory with inconsistent t and Δt initially
     N = 10
     Δt_val = 0.5
 
     traj = NamedTrajectory(
         (
-            x = rand(2, N),
-            u = rand(1, N),
+            x = rand(rng, 2, N),
+            u = rand(rng, 1, N),
             Δt = fill(Δt_val, N),
-            t = cumsum(rand(N)),  # Random times - inconsistent!
+            t = cumsum(rand(rng, N)),  # Random times - inconsistent!
         );
         controls = (:u, :Δt, :t),  # t is also a control to be optimized
         timestep = :Δt,
@@ -123,3 +129,73 @@ end
     # Verify initial time constraint
     @test abs(t[1]) < 1e-8
 end
+
+@testitem "TimeConsistencyConstraint free-time robustness sweep" begin
+    include("../../../test/test_utils.jl")
+    using NamedTrajectories
+    using Random
+
+    # K independent seeds; pass if ≥80% land within the same tolerance the
+    # deterministic test uses. Catches regressions where the solver's local-
+    # minimum behavior degrades for "typical" inconsistent initializations.
+    function run_sweep(K::Int)
+        pass_count = 0
+        failures = String[]
+        for seed = 1:K
+            rng = MersenneTwister(seed)
+            N = 10
+            Δt_val = 0.5
+
+            traj = NamedTrajectory(
+                (
+                    x = rand(rng, 2, N),
+                    u = rand(rng, 1, N),
+                    Δt = fill(Δt_val, N),
+                    t = cumsum(rand(rng, N)),
+                );
+                controls = (:u, :Δt, :t),
+                timestep = :Δt,
+                bounds = (u = (-1.0, 1.0), t = (0.0, 10.0)),
+                initial = (t = [0.0],),
+            )
+
+            J = QuadraticRegularizer(:u, traj, 1.0)
+            J += QuadraticRegularizer(:t, traj, 0.1)
+
+            time_con = TimeConsistencyConstraint()
+
+            prob = DirectTrajOptProblem(
+                traj,
+                J,
+                AbstractIntegrator[];
+                constraints = [time_con],
+            )
+
+            ok = try
+                solve!(prob; max_iter = 100)
+                t = prob.trajectory.t
+                Δt = prob.trajectory.Δt
+                consistent = all(k -> abs(t[k+1] - t[k] - Δt[k]) < 1e-6, 1:(N-1))
+                initial_ok = abs(t[1]) < 1e-8
+                consistent && initial_ok
+            catch e
+                push!(failures, "seed $seed: threw $(typeof(e))")
+                false
+            end
+
+            if ok
+                pass_count += 1
+            else
+                push!(failures, "seed $seed: tolerance not met after solve")
+            end
+        end
+        return pass_count, failures
+    end
+
+    K = 20
+    pass_threshold = 0.80
+    pass_count, failures = run_sweep(K)
+    pass_rate = pass_count / K
+    @info "TimeConsistencyConstraint free-time robustness sweep" pass_count K pass_rate failures
+    @test pass_rate >= pass_threshold
+end

src/constraints/nonlinear/knot_point_constraint.jl

@@ -312,12 +312,17 @@ end
 
 @testitem "NonlinearKnotPointConstraint - single variable with vector syntax" begin
     using DirectTrajOpt: CommonInterface
+    using Random
 
     include("../../../test/test_utils.jl")
 
-    _, traj = bilinear_dynamics_and_trajectory()
+    # Deterministic trajectory + multiplier: the test verifies vector vs
+    # single-symbol syntax structurally, plus a single finite-diff comparison.
+    # Multi-seed coverage of the finite-diff tolerance lives in the robustness
+    # testitem below — this one must be reproducible across Julia versions.
+    rng = MersenneTwister(0)
+    _, traj = bilinear_dynamics_and_trajectory(; rng = rng)
 
-    # Test that [:u] syntax works the same as :u
     g(a) = [norm(a) - 1.0]
 
     NLC1 = NonlinearKnotPointConstraint(g, :u, traj; equality = false)
@@ -330,9 +335,78 @@ end
 
     @test δ1 ≈ δ2
 
-    # Test both with finite differences
-    test_constraint(NLC1, traj; atol = 1e-3)
-    test_constraint(NLC2, traj; atol = 1e-3)
+    # Test both with finite differences (seeded μ → deterministic Hessian check)
+    test_constraint(NLC1, traj; atol = 1e-3, rng = MersenneTwister(1))
+    test_constraint(NLC2, traj; atol = 1e-3, rng = MersenneTwister(1))
+end
+
+@testitem "NonlinearKnotPointConstraint vector syntax robustness sweep" begin
+    using DirectTrajOpt: CommonInterface
+    using Random
+
+    include("../../../test/test_utils.jl")
+
+    # Finite-diff vs analytic Jacobian/Hessian at atol=1e-3 is fundamentally
+    # noisy in the random trajectory point. Run K seeds; require ≥80% pass.
+    # `assert=false`: inspect jacobian_pass/hessian_pass directly instead of
+    # registering @test outcomes per-seed. `test_equality=false` uses the
+    # norm-based aggregate check (single noisy entries shouldn't fail the
+    # seed — the overall derivative match is what matters).
+    function run_sweep(K::Int)
+        pass_count = 0
+        failures = String[]
+        g(a) = [norm(a) - 1.0]
+        for seed = 1:K
+            _, traj = bilinear_dynamics_and_trajectory(; rng = MersenneTwister(seed))
+            NLC1 = NonlinearKnotPointConstraint(g, :u, traj; equality = false)
+            NLC2 = NonlinearKnotPointConstraint(g, [:u], traj; equality = false)
+
+            δ1 = zeros(NLC1.dim)
+            δ2 = zeros(NLC2.dim)
+            CommonInterface.evaluate!(δ1, NLC1, traj)
+            CommonInterface.evaluate!(δ2, NLC2, traj)
+            syntax_ok = isapprox(δ1, δ2)
+
+            r1 = test_constraint(
+                NLC1,
+                traj;
+                atol = 1e-3,
+                test_equality = false,
+                rng = MersenneTwister(seed + K),
+                assert = false,
+            )
+            r2 = test_constraint(
+                NLC2,
+                traj;
+                atol = 1e-3,
+                test_equality = false,
+                rng = MersenneTwister(seed + K),
+                assert = false,
+            )
+            fd_ok =
+                r1.jacobian_pass && r1.hessian_pass && r2.jacobian_pass && r2.hessian_pass
+
+            if syntax_ok && fd_ok
+                pass_count += 1
+            else
+                !syntax_ok && push!(failures, "seed $seed: δ1 ≉ δ2 (syntax check)")
+                !fd_ok && push!(failures, "seed $seed: finite-diff check failed")
+            end
+        end
+        return pass_count, failures
+    end
+
+    # FD-vs-analytic Jacobian/Hessian on a norm comparison is inherently
+    # noisier than the time-consistency check. Observed K=20 pass rate on
+    # Julia 1.12 is 0.80; threshold relaxed to 0.65 so a one-seed shift
+    # across Julia versions doesn't false-fail the gate. A true regression
+    # in the analytic derivative would drop the rate well below 0.5.
+    K = 20
+    pass_threshold = 0.65
+    pass_count, failures = run_sweep(K)
+    pass_rate = pass_count / K
+    @info "NonlinearKnotPointConstraint vector-syntax robustness sweep" pass_count K pass_rate failures
+    @test pass_rate >= pass_threshold
 end
 
 @testitem "NonlinearKnotPointConstraint - multiple variables concatenated" begin

test/runtests.jl

@@ -3,5 +3,7 @@ using TestItemRunner
 
 include("test_snippets.jl")
 
-# Run all testitem tests in package
-@run_package_tests
+# Exclude `benchmark/` testitems from the main test run — they live in a
+# subproject with its own Project.toml (different deps, e.g. HarmoniqsBenchmarks)
+# and are exercised by a dedicated workflow.
+@run_package_tests filter = ti -> !occursin("/benchmark/", ti.filename)