From dbcc207f7b0bf7f2a44c2ccad8fd0113aef6c4fd Mon Sep 17 00:00:00 2001 From: abavoil <17646791+abavoil@users.noreply.github.com> Date: Thu, 3 Sep 2026 16:25:05 +0200 Subject: [PATCH] src/bvp: fix uninitialized residual tail and size contract inconsistencies bvp_residual for Shooting/Trapeze/Collocation allocated out = similar(X) but only wrote the prefix out[1:n*M]; the trailing 'period' slot (advertised by generate_solution and total_dim) was never written, so residual(prob, x, p) returned uninitialized memory in newton (nondeterministic results) and the AutoDiffDense Jacobian carried a garbage last row. Fixed-interval design, consistent with #312 ('no need for phase') and #315 ('fixed final time; any unknown like tf must be a constant state'): - residuals fully write an output of length solution_dim and reject wrong-size input with an ArgumentError - generate_solution returns n*M (no dead trailing slot) - Base.length(DiscretizedBVP) = solution_dim; total_dim removed - get_periodic_orbit uses the model time interval instead of X[end] - remove dead bvp_jacobian(::Nothing) for Shooting (called three undefined helpers: _shooting_phase, integrate_shooting, euler_integrate) and the dangling exports - tests: update test/bvp/bvp.jl to the solution_dim contract, add dedicated MRE regression test test/bvp/mre_shooting_uninit.jl --- src/bvp/BVP.jl | 4 -- src/bvp/DiscretizedBVP.jl | 4 +- src/bvp/Discretizers.jl | 5 +- src/bvp/collocation/residual.jl | 12 ++-- src/bvp/discretize.jl | 8 +-- src/bvp/integration.jl | 6 +- src/bvp/shooting/jacobian.jl | 105 -------------------------------- src/bvp/shooting/residual.jl | 15 +++-- src/bvp/trapeze/residual.jl | 11 ++-- test/bvp/bvp.jl | 50 ++++++++++++--- test/bvp/mre_shooting_uninit.jl | 49 +++++++++++++++ 11 files changed, 124 insertions(+), 145 deletions(-) delete mode 100644 src/bvp/shooting/jacobian.jl create mode 100644 test/bvp/mre_shooting_uninit.jl diff --git a/src/bvp/BVP.jl b/src/bvp/BVP.jl index d5e2f934..d0c11d46 100644 --- a/src/bvp/BVP.jl +++ b/src/bvp/BVP.jl @@ -59,7 +59,6 @@ include("discretize.jl") # Residual/Jacobian implementations for each discretizer include("shooting/residual.jl") -include("shooting/jacobian.jl") # Shooting has specialized analytical jacobian include("trapeze/residual.jl") include("trapeze/jacobian.jl") include("collocation/residual.jl") @@ -83,7 +82,4 @@ export state_dimension, getperiod export BVPBifProblem export get_periodic_orbit, get_bvp -# Internal exports for extensions -export integrate_shooting, integrate_with_sensitivity - end # module BVP diff --git a/src/bvp/DiscretizedBVP.jl b/src/bvp/DiscretizedBVP.jl index 49afd409..762e73c8 100644 --- a/src/bvp/DiscretizedBVP.jl +++ b/src/bvp/DiscretizedBVP.jl @@ -78,8 +78,8 @@ end """State dimension.""" state_dimension(bvp::DiscretizedBVP) = state_dimension(bvp.model) -"""Total dimension of the discretized problem.""" -Base.length(bvp::DiscretizedBVP) = total_dim(bvp.discretizer, state_dimension(bvp)) +"""Total dimension of the discretized problem (the time span is fixed by the model, there is no trailing period slot).""" +Base.length(bvp::DiscretizedBVP) = solution_dim(bvp.discretizer, state_dimension(bvp)) """Get the underlying model.""" get_model(bvp::DiscretizedBVP) = bvp.model diff --git a/src/bvp/Discretizers.jl b/src/bvp/Discretizers.jl index e47b20ca..9cd70197 100644 --- a/src/bvp/Discretizers.jl +++ b/src/bvp/Discretizers.jl @@ -185,15 +185,12 @@ mesh_size(d::Shooting) = d.M mesh_size(d::Trapeze) = d.M mesh_size(d::Collocation) = d.Ntst * d.m + 1 -"""Dimension of the discretized solution vector (excluding period).""" +"""Dimension of the discretized solution vector.""" function solution_dim end solution_dim(d::Shooting, n::Int) = n * d.M solution_dim(d::Trapeze, n::Int) = n * d.M solution_dim(d::Collocation, n::Int) = n * (d.Ntst * d.m + 1) - -"""Total dimension including period/parameter.""" -total_dim(d::AbstractDiscretizer, n::Int) = solution_dim(d, n) + 1 #━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ # Display #━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ diff --git a/src/bvp/collocation/residual.jl b/src/bvp/collocation/residual.jl index 38a62251..b3aeeb7f 100644 --- a/src/bvp/collocation/residual.jl +++ b/src/bvp/collocation/residual.jl @@ -8,13 +8,15 @@ function bvp_residual(d_bvp::DiscretizedBVP{<: BVPModel, <: Collocation}, X, p) interval = get_time_interval(model) δT = interval[2] - interval[1] + N = nf * N_total + length(X) == N || throw(ArgumentError("bvp_residual: expected length(X) == $N, got $(length(X))")) + # Extract solution - Xm = reshape(@view(X[1:nf*N_total]), nf, N_total) + Xm = reshape(@view(X[1:N]), nf, N_total) - # Get output buffer from cache - # Robust check: only use cache for Float64 to avoid chunk mismatch in Dual - out = similar(X) - outm = reshape(@view(out[1:nf*N_total]), nf, N_total) + # Allocate output; every entry is written below + out = similar(X, N) + outm = reshape(out, nf, N_total) # Core residual computation from BifurcationKit # This writes to outm[:, 1:Ntst*m] diff --git a/src/bvp/discretize.jl b/src/bvp/discretize.jl index b3f1e3e8..801d9833 100644 --- a/src/bvp/discretize.jl +++ b/src/bvp/discretize.jl @@ -150,8 +150,8 @@ function generate_solution(model::BVPModel, disc::Shooting, cache, orbit) M = disc.M t0, tf = get_time_interval(model) T = tf - t0 - # Sample at M shooting points - X = zeros(n * M + 1) + # Sample at M shooting points; the time span is fixed by the model + X = zeros(n * M) for i in 1:M t = (i - 1) / M * T X[(i-1)*n+1 : i*n] .= orbit(t) @@ -164,8 +164,8 @@ function generate_solution(model::BVPModel, disc::Trapeze, cache, orbit) M = disc.M t0, tf = get_time_interval(model) T = tf - t0 - # Sample at M time slices - X = zeros(n * M + 1) + # Sample at M time slices; the time span is fixed by the model + X = zeros(n * M) for i in 1:M t = (i - 1) / (M - 1) * T X[(i-1)*n+1 : i*n] .= orbit(t) diff --git a/src/bvp/integration.jl b/src/bvp/integration.jl index 1b9f76f1..322d6bdf 100644 --- a/src/bvp/integration.jl +++ b/src/bvp/integration.jl @@ -22,8 +22,10 @@ Returns a `NamedTuple` with fields: function get_periodic_orbit(bvp::DiscretizedBVP, X, p) n = state_dimension(bvp) disc = get_discretizer(bvp) - T = X[end] - + # The time span is fixed by the model, not stored in X + t0, tf = get_time_interval(get_model(bvp)) + T = tf - t0 + return _get_periodic_orbit(disc, X, n, T) end diff --git a/src/bvp/shooting/jacobian.jl b/src/bvp/shooting/jacobian.jl deleted file mode 100644 index 2d48bc12..00000000 --- a/src/bvp/shooting/jacobian.jl +++ /dev/null @@ -1,105 +0,0 @@ -# Shooting Jacobian Implementation - -import BifurcationKit: AbstractJacobianType, AutoDiffDense - -""" -$(TYPEDSIGNATURES) - -Compute the Jacobian for shooting discretization. - -## Jacobian Structure (for M=3 intervals, periodic): - ┌ ┐ - │ Φ₁ -I 0 ∂R₁/∂T │ - │ 0 Φ₂ -I ∂R₂/∂T │ - │ -I 0 Φ₃ ∂R₃/∂T │ - │ ∂g/∂u₁ 0 0 ∂g/∂T │ - └ ┘ - -Where Φᵢ = ∂φ(uᵢ)/∂uᵢ is the monodromy matrix. -""" - -function bvp_jacobian(bvp::DiscretizedBVP{<:BVPModel, <:Shooting}, jac::Nothing, X, p) - model = bvp.model - disc = bvp.discretizer - - n = state_dimension(model) - M = disc.M - N = n * M + 1 - - # Extract shooting points and period - U = reshape(@view(X[1:n*M]), n, M) - T = X[end] - dt = T / M - - # Allocate Jacobian - J = zeros(eltype(X), N, N) - In = Matrix{eltype(X)}(LinearAlgebra.I, n, n) - - for i in 1:M - uᵢ = U[:, i] - row_start = (i - 1) * n + 1 - row_end = i * n - - # Compute flow and monodromy matrix - u_final, Φᵢ = integrate_with_sensitivity(model.F, uᵢ, p, dt, disc.alg) - - # ∂Rᵢ/∂uᵢ = Φᵢ (monodromy matrix) - col_start = (i - 1) * n + 1 - col_end = i * n - J[row_start:row_end, col_start:col_end] .= Φᵢ - - # ∂Rᵢ/∂uᵢ₊₁ = -I - i_next = (i == M) ? 1 : i + 1 - col_start_next = (i_next - 1) * n + 1 - col_end_next = i_next * n - J[row_start:row_end, col_start_next:col_end_next] .= -In - - # ∂Rᵢ/∂T = F(u_final) / M - J[row_start:row_end, end] .= model.F(u_final, p) ./ M - end - - # Last row: phase condition Jacobian - u0 = U[:, 1] - uM = U[:, M] - uT = integrate_shooting(model.F, uM, p, dt, disc.alg) - - # ∂g/∂u₁ - J[end, 1:n] .= ForwardDiff.gradient(u -> _shooting_phase(model, u, uT, p), u0) - - # ∂g/∂T (finite difference for simplicity) - δ = 1e-8 - uT_plus = integrate_shooting(model.F, uM, p, dt + δ/M, disc.alg) - phase_plus = _shooting_phase(model, u0, uT_plus, p) - phase_current = _shooting_phase(model, u0, uT, p) - J[end, end] = (phase_plus - phase_current) / δ - - return J -end - -""" -Integrate ODE and return both final state and monodromy matrix. - -This is a dispatch point: -- `alg=nothing`: Uses Euler + ForwardDiff for monodromy -- `alg::OrdinaryDiffEqAlgorithm`: Uses variational equations (requires extension) -""" -function integrate_with_sensitivity end - -# Default: nothing algorithm uses Euler + AD -function integrate_with_sensitivity(F, u0, p, dt, ::Nothing) - # Euler integration + AD for monodromy matrix - u_final = euler_integrate(F, u0, p, dt) - Φ = ForwardDiff.jacobian(u -> euler_integrate(F, u, p, dt), u0) - return u_final, Φ -end - -# Fallback for unknown algorithm types -function integrate_with_sensitivity(F, u0, p, dt, alg) - error(""" - Sensitivity computation with algorithm $(typeof(alg)) requires OrdinaryDiffEq.jl. - - Either: - 1. Use `Shooting(alg=nothing)` for built-in Euler integration - 2. Load OrdinaryDiffEq: `using OrdinaryDiffEq` - """) -end diff --git a/src/bvp/shooting/residual.jl b/src/bvp/shooting/residual.jl index 4109a697..cbe347cd 100644 --- a/src/bvp/shooting/residual.jl +++ b/src/bvp/shooting/residual.jl @@ -11,14 +11,19 @@ function bvp_residual(d_bvp::DiscretizedBVP{<:BVPModel, <:Shooting}, X, p) t0, tf = get_time_interval(model) M = mesh_size(disc) + # The time span is fixed by the model (free quantities like tf must be + # encoded as constant state components, cf. issues #312/#315) + N = n * M + length(X) == N || throw(ArgumentError("bvp_residual: expected length(X) == $N, got $(length(X))")) + # Extract shooting points and period - Xm = reshape(@view(X[1:n*M]), n, M) + Xm = reshape(@view(X[1:N]), n, M) T = tf - t0 - # Allocate output - out = similar(X) - outm = reshape(@view(out[1:n*M]), n, M) - + # Allocate output; every entry is written by bvp_residual_bare! + out = similar(X, N) + outm = reshape(out, n, M) + # Core residual computation using BVP-specific po_residual_bare! bvp_residual_bare!(d_bvp, outm, Xm, p, T) return out diff --git a/src/bvp/trapeze/residual.jl b/src/bvp/trapeze/residual.jl index 5c70f0e9..024f4477 100644 --- a/src/bvp/trapeze/residual.jl +++ b/src/bvp/trapeze/residual.jl @@ -5,15 +5,18 @@ function bvp_residual(d_bvp::DiscretizedBVP{<: BVPModel, <: Trapeze}, X, p) n = state_dimension(model) M = disc.M + N = n * M + length(X) == N || throw(ArgumentError("bvp_residual: expected length(X) == $N, got $(length(X))")) + interval = get_time_interval(model) δT = interval[2] - interval[1] # Extract time slices - Xm = reshape(@view(X[1:n*M]), n, M) + Xm = reshape(@view(X[1:N]), n, M) - # Get output buffer; element type follows X (covers Dual numbers in AD) - out = similar(X) - outm = reshape(@view(out[1:n*M]), n, M) + # Allocate output; every entry is written below (element type follows X, covers Dual numbers in AD) + out = similar(X, N) + outm = reshape(out, n, M) # Sequential trapezoid scheme for M-1 intervals using potrap_scheme!. # M time points span [t0, tf] via M-1 intervals with normalized mesh weights. diff --git a/test/bvp/bvp.jl b/test/bvp/bvp.jl index 0870c99e..f089c751 100644 --- a/test/bvp/bvp.jl +++ b/test/bvp/bvp.jl @@ -1,4 +1,5 @@ using BifurcationKit, Test +import OrdinaryDiffEq as ODE const BK = BifurcationKit # ============================================================================== @@ -56,7 +57,6 @@ let @test disc.M == 4 @test BK.BVP.mesh_size(disc) == 4 @test BK.BVP.solution_dim(disc, 2) == 8 - @test BK.BVP.total_dim(disc, 2) == 9 end # ----- Trapeze ----- @@ -66,7 +66,6 @@ let @test disc.M == 50 @test BK.BVP.mesh_size(disc) == 50 @test BK.BVP.solution_dim(disc, 2) == 100 - @test BK.BVP.total_dim(disc, 2) == 101 mesh = BK.TimeMesh([0.5, 0.5]) disc2 = BK.BVP.Trapeze(; M=3, mesh) @@ -83,7 +82,6 @@ let @test disc.m == 4 @test BK.BVP.mesh_size(disc) == 20*4 + 1 @test BK.BVP.solution_dim(disc, 2) == 2 * (20*4 + 1) - @test BK.BVP.total_dim(disc, 2) == 2 * (20*4 + 1) + 1 end # ----- discretize: Trapeze ----- @@ -99,7 +97,7 @@ let @test bvp.model === model @test bvp.discretizer === disc @test BK.BVP.state_dimension(bvp) == 2 - @test length(bvp) == 2*10 + 1 + @test length(bvp) == 2*10 end # ----- discretize: Collocation ----- @@ -112,7 +110,7 @@ let show(bvp) @test bvp isa BK.BVP.DiscretizedBVP @test BK.BVP.state_dimension(bvp) == 2 - @test length(bvp) == 2 * (5*3 + 1) + 1 + @test length(bvp) == 2 * (5*3 + 1) end # ----- DiscretizedBVP getters ----- @@ -170,7 +168,7 @@ let x0 = zeros(length(bvp)) res = BK.BVP.bvp_residual(bvp, x0, (ω=1.0,)) @test length(res) == length(bvp) - @test res[1:end-1] == zeros(length(bvp)-1) + @test res == zeros(length(bvp)) # fully written, no uninitialized tail end # ----- bvp_jacobian: Trapeze, FD vs analytical ----- @@ -303,9 +301,10 @@ let bvp = BK.BVP.discretize(model, disc) x0 = BK.BVP.generate_solution(bvp, t -> [cos(t), sin(t)]) @test length(x0) == length(bvp) - @test x0[end] == 0.0 # T = 0 in initial guess # First slice should be orbit at t=0 @test x0[1:2] ≈ [1.0, 0.0] + # Last slice is orbit(T) with T = tf - t0 = 1 + @test x0[end] == sin(1.0) end # ----- generate_solution (Collocation) ----- @@ -452,7 +451,7 @@ let @test haskey(po, :period) @test length(po.t) == disc.M @test size(po.u, 2) == disc.M - @test po.period == x0[end] + @test po.period == 1.0 # fixed time interval (0, 1) of the model end # ----- get_periodic_orbit (Collocation) ----- @@ -469,7 +468,7 @@ let @test haskey(po, :u) @test haskey(po, :period) @test length(po.t) == disc.Ntst * disc.m + 1 - @test po.period == x0[end] + @test po.period == 1.0 # fixed time interval (0, 1) of the model end # ----- update_phase_reference! ----- @@ -482,7 +481,7 @@ let @test BK.BVP.update_phase_reference!(bvp, rand(length(bvp)), (ω=1.0,)) == true end -# ----- Collocation: bvp_jacobian with AutoDiffDense (default, expects full vector with period) ----- +# ----- Collocation: bvp_jacobian with AutoDiffDense (default) ----- let F(u, p) = [u[2], -p.ω^2 * u[1]] g(u0, uT, p) = [u0[1], uT[1]] @@ -521,3 +520,34 @@ let @test br isa BK.AbstractBranchResult @test length(br) > 0 end + +# ----- Regression: Shooting residual is fully written, Newton is deterministic ----- +# (the residual used to leave its last entry uninitialized: `out = similar(X)` with +# only the prefix out[1:n*M] written; the trailing +1 "period" slot was never filled) +let + F(u, p, t=0) = [u[2], -p.μ * u[1]] + g(u0, uT, p) = [u0[1] - 1.0, u0[2]] # pins the phase, the system is square + n, M = 2, 4 + + odeprob = ODE.ODEProblem(F, [1.0, 0.0], (0.0, 2π), (μ = 1.0,)) + model = BK.BVP.BVPModel(odeprob, g; n) + bvp = BK.BVP.discretize(model, BK.BVP.Shooting(M, ODE.Tsit5(), false)) + + x0 = BK.BVP.generate_solution(bvp, t -> [cos(t), sin(t)]) + @test length(x0) == n * M + x0[1] += 0.05 # perturb so that Newton actually iterates + + prob = BK.BVP.BVPBifProblem(bvp, x0, (μ = 1.0,), (@optic _.μ)) + + # identical inputs must give identical residuals (no uninitialized memory) + r1 = BK.residual(prob, prob.u0, prob.params) + r2 = BK.residual(prob, prob.u0, prob.params) + @test r1 == r2 + @test length(r1) == n * M + + # wrong-size input must fail loudly instead of reading garbage + @test_throws ArgumentError BK.BVP.bvp_residual(bvp, vcat(x0, 0.0), (μ = 1.0,)) + + sol = BK.solve(prob, BK.Newton(), NewtonPar(tol = 1e-9, verbose = false)) + @test BK.converged(sol) +end diff --git a/test/bvp/mre_shooting_uninit.jl b/test/bvp/mre_shooting_uninit.jl new file mode 100644 index 00000000..ce393959 --- /dev/null +++ b/test/bvp/mre_shooting_uninit.jl @@ -0,0 +1,49 @@ +using BifurcationKit, Test +import OrdinaryDiffEq as ODE +const BK = BifurcationKit + +# ============================================================================== +# Regression test for the uninitialized residual tail of the BVP Shooting +# discretization. +# +# bvp_residual used to allocate `out = similar(X)` of length n*M+1 (trailing +# "period" slot advertised by generate_solution/total_dim) but wrote only +# out[1:n*M], so `residual(prob, x, p)` — read at Newton.jl:74 — returned +# uninitialized memory (NaN / subnormals, differing between identical calls) +# and the AutoDiffDense Jacobian carried a garbage last row. Newton was +# nondeterministic. +# +# Fix (fixed-interval design): the time span is fixed by the BVPModel, there is +# no trailing slot; residuals fully write an output of length solution_dim and +# reject wrong-size input. Unknowns like tf must be constant state components. +# ============================================================================== + +let + # harmonic oscillator on [0, 2π], periodic boundary conditions + F(u, p, t=0) = [u[2], -p.μ * u[1]] + g(u0, uT, p) = u0 .- uT + n, M = 2, 4 + + odeprob = ODE.ODEProblem(F, [1.0, 0.0], (0.0, 2π), (μ = 1.0,)) + model = BK.BVP.BVPModel(odeprob, g; n) + bvp = BK.BVP.discretize(model, BK.BVP.Shooting(M, ODE.Tsit5(), false)) + + x0 = BK.BVP.generate_solution(bvp, t -> [cos(t), sin(t)]) + @test length(x0) == n * M # no trailing period slot + + prob = BK.BVP.BVPBifProblem(bvp, x0, (μ = 1.0,), (@optic _.μ)) + + # residual must be fully written: identical inputs, identical outputs + r1 = BK.residual(prob, prob.u0, prob.params) + r2 = BK.residual(prob, prob.u0, prob.params) + @test length(r1) == n * M + @test r1 == r2 + + # wrong-size input must fail loudly instead of reading uninitialized memory + @test_throws ArgumentError BK.BVP.bvp_residual(bvp, vcat(x0, 0.0), (μ = 1.0,)) + @test_throws ArgumentError BK.BVP.bvp_residual(bvp, x0[1:end-1], (μ = 1.0,)) + + # Newton on the Shooting path is deterministic and converges + sol = BK.solve(prob, BK.Newton(), NewtonPar(tol = 1e-10, verbose = false)) + @test BK.converged(sol) +end