src/bvp: fix uninitialized residual tail and size contract inconsistencies

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

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

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
 #━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

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]

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)

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
 

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

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.

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